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Part I 
Lectures 



Lecture 1 

Basic Quantum Mechanics 



1.1 Operators in quantum mechanics 

An observable is a dynamical, measurable variable of a system. In classical 
physics observables are described using functions, in quantum mechanics 
operators are used. An operator is a symbol for a recipe how an operation 
on a function has to be performed. 

1.1.1 Eigenfunctions and eigenvalues 

In general an operator creates a new function when it is applied to a function. 
In special cases the new function is the same as the original one multiplied 
by a constant. A function / is an eigenfunction of an operator, if 

Clf = uf , 

where u is a constant. Such an equation is called an eigenvalue equation, u 
is the eigenvalue of the operator Cl. 

The set of all eigenfunctions of an operator is a complete basis, this means an 
arbitrary function can be written as a linear combination of eigenfunctions. 

A function g is expanded in eigenfunctions 

n 

with c n the expansion coefficients and the sum runs over all functions /„. We 
can easily calculate the effect of the operator Cl on the function g 

Clg Cl ^ ^ Cnfn ^ ^ C n Clf n ^ ^ C n U n f n 
n n n 
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Quantum mechanics can be developed 
However, it is instructive and intuitive 
most common uses the following form 
(position representation) 

The special form of the momentum 
introduced in the next section. 



without explicitely defining operators, 
to use special forms of operators. The 
for position and momentum operator 

h_d_ 

^ x 2ni dx 

operator follows from the postulates 



1.1.2 Commutators 

The order in which operators are applied to a function is important. Gener- 
ally, we have 

ABf ^ BAf , 

or simply AB ^ BA. The quantity AB — BA is called the commutator of A 
and B. 

[A, B]= AB-BA . 
Operators with the same set of eigenfunctions commute 

ABg = BAg 

A^UJ^Cnfn = B^UJ^Cnfn 

n n 

EAR r \ v R R r 

U n UJ n C n f n = 2_^U) n UJ n C n f n 

n n 

The opposite is also true; operators that commute can be diagonalized with 
the same set of functions. 

The commutator for position and momentum operator in the special repre- 
sentation introduced above is 

r , v, h df h d 

\X, Vx \t = (XPx — PxX) t — X ■ — — — • XT 

l ,yx\j v fx fx j j 2m dx 2m dx 
2ni dx 2-ni 2ni dx 

ih 
~ 2tt '"' 

This is true for any function /, therefore we write 

r i ih 
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1.1.3 Construction of operators 

Many operators we will need can be constructed from the operators for 
position and momentum. For example the operator of the kinetic energy 
T = p 2 /2m is in one dimension 

f= vl i fhd\ 2 h 2 d 2 



2m 2m \ i dx J 2m dx 2 

and in three dimensions 



T 



Vl +Pl+P 2 z 

2m 

a 2 r d 2 



+ 







2m 1 


\l ox 


d 2 
dy 2 


d 2 1 


dz 2 J 




2m \ dx 2 
2m 

where h — h/2n and the operator V 2 is the Laplacian. 
The potential energy of a particle V(x) is a multiplicative function in the 
position representation. This is also the case for the operator of the potential 
energy, for example the operator of the Coulomb energy of an electron in the 
field of a nuclei with charge Z 

v = _ *L. 

The operator of the total energy of a system is the Hamilton operator 

H = f + V 

The general recipe for the construction of quantum mechanical operators is: 

1. Write down the classical expression for the observable using the position 
and momentum variables of the system. 

2. Replace the different functions with their quantum mechanical opera- 
tors. 



1.2 Postulates of quantum mechanics 
1.2.1 Postulates 

The following postulates are not complete and in their most strict form, but 
they are sufficient for our purposes. 
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Table 1.1: Transition from classical expressions to quantum mechanical op- 
erators. 



classical 



quantum mechanical 



Position 
e.g. x, y 

Functions of position 
e.g. x 2 y 

or general f(x,y,z) 

Momentum 
e.g. p x 

Functions of momentum 
e.g. (p*) 2 

or general g(p x ,p y ,p z ) 

Functions of position and momentum 
e.g. I = r x p 

Time t 

Energy E 



x- , y- 

multiplicative operators 

2 

x y- 

f(x,y,z)- 



hd_ 

i dx 

differential operators 

v i dx ' i dyi i dz> 



(x,y,z) X | (j^, ^, 
t- 

UL at 



Postulate 1 

The state of a system is completely characterized by its wave function 

*(n,r 2 ,...,*)- 

The positions r±, r 2 , . . . are associated with particles 1, 2, ... . In general 
additional internal coordinates (spin) are needed. From the wave function 
all information on the properties of a system can be derived, that are also 
accessible by experiments. 

Postulate 2 

Observables are represented by operators which have the following prop- 
erties 

[q,Pq>] = ih5 m , , [q,q] = , \p q ,P q '\ = , 

where q and q' are one of the coordinates x, y, z and p q and p q i are the 
corresponding momenta. 
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With this postulate all operators that are functions of position and momen- 
tum coordinates are defined. Internal (spin) coordinates have to be defined 
separately. 

Postulate 3 

For a system characterized by its wave function \l/ the average value of 
a series of measurements of the observable Q is given by the expectation 
value of the corresponding operator. 

The expectation value of an operator Cl for a general state \I/ is defined by 
If the wave function is normalized 



"Vdr = 1 , 

we get for the expectation value 

(Q) = J V*£wdr . 

We will assume that state functions are normalized. If \I/ is an eigenstate of 
Q with eigenvalue uj, we have 

(Q) = J y*Q^dT = J m*uumdr = uo J m*mdr = uo . 

A series of measurements on identical systems will give the averaged value 
uj. Assuming that \1/ is an eigenstate of the system, but not of Cl, we can 
expand \1/ in the eigenfunctions of Cl 

^ = c n^n , where tt^ n = uj n ^ n gilt. 

n 

The expectation value is 

(ft) = [Yl Cm ^ m J n [Yl Cn ^ n ) dr = Yl C ™° n / ^m^ndT 

\ m J \ n J nm 



J2c* m c n u n J ^ n dr 



Because the eigenfunctions are an orthonormal basis, the last integral van- 
ishes for m 7^ n and we get 

(^) = ^2 C m C n^n J ^* rn ^ ' n^T = ^ C* n C n U n = ^\c n \ 2 UJ n . 

nm n n 

The expectation value is therefore a weighted sum of the eigenvalues of Q. 
The contribution of a certain eigenvalue is given by the square of the absolute 
value of the corresponding expansion coefficient. 

Postulate 3' 

If \& is an eigenfunction of Cl all measurements will give the value uj. Is 
\1> not an eigenfunction of Cl one will get with each measurement a value 
u n , an eigenvalue of Cl. The probability to get a certain eigenvalue is 
proportional to |c„| 2 , where c„ is the overlap of \I/ with the eigenfunction 



\& n of Q. 



Postulate 4 

The probability to find a particle in the volume element dr at position r 
is proportional to \^/(r)\ 2 dr. 

This probability interpretation of the wave function is due to Max Born. 
From this it follows that ^(r)! 2 is a probability density. The wave function 
is a probability amplitude. From this interpretation it also follows that the 
wave functions have to be square integrable 



\% 2 dr < oo , 

Because the probability to find a particle at any place in space has to be 
finite. 

Postulate 5 

The time evolution of a wave function \&(ri, r 2 , . . . , t) is governed by the 
Schrodinger equation 

ih— = Hm . 

dt 

This is a partial differential equation. H is the Hamilton operator of the 
system. For example the Schrodinger equation for a particle in one dimension 
is 

H 2 d 2 ^ _ _ / NT 
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1.2.2 Separation of the Schrodinger equation 

Often the Schrodinger equation can be separated into an equation for the 
positions and an equation for the time evolution. If the Hamilton operator 
does not depend on time, we get in one dimension 

m = "^ a? + y W* = B * = m m 

For stationary states the energy is a constant and we can solve the ordinary 
differential equation using 

(E 

^(x, t) = ip(x) exp —t 
Inserting this into the Schrodinger equation we get 

the time-independent Schrodinger equation which determines the function 

ip(x). 



1.3 Hydrogen Atom 

The hydrogen atom consists of a proton of mass m p and an electron of mass 
m e . Reducing the two particle problem to a single particle equation can be 
achieved by separation of the center of mass movement: 

• 3 center of mass coordinates (translation) 

• 3 inner coordinates (relative movement) 

This separation leads to the reduced mass for the relative movement: 

fx = e —^- = 0.999456m e w m e 

m e + m p 

H = -^V 2 + V(r) 

The relative movement can be cast into spherical coordinates (r, 0, 0). The 
Coulomb potential for this special case is a central field problem 

V(r)= 



r 

The nuclear charge is Z = 1 for the hydrogen atom. 
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1.3.1 Angular momentum operators 

1.3.2 Definitions 

Classical angular momentum: 

I = j- x p = (l x , l y , l z ) 
l x = VPz ~ zp y 

ly Xp z 

l z = xp y - yp x 
Operators in quantum mechanics: 

f h ( d d 

lx = -\ Vtt - z— 

% \ oz oy 

f _ h f d_ _ d_ 



I _ h ( d _ d 
i \ dy ^ dx 



Total angular momentum : 
classical: 



1^=1.1 = % + % + % 



in quantum mechanics: 



I l x -\- ly -\- l z Ixl'x lyly Iz^z 



Commutators 

For two components: 

iJ y f = - h2 { y w z -%)( z Tx- x h f 



h 2 , d 2 f M df d 2 f 2 d 2 f d 2 f 
-h yz—— — h y— yx-—r — z — — — h zx- 



dzdx dx dz 2 dydx dydz ) 

. , ■ d 2 f df d 2 f 2 iff iff \ 

lylxf = -n [zy^- + x— -xy— - z —— + xz—— ) 
dxdz ay oz 1 dxdy ozoy J 



Differences 



(ij y - l y l x )f = -h 2 {y^ - 



h ( df df 

= ih- x- y— 

% \ ay ox 

= itU e f 
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Commutator-equations: 



For a single component and I 2 : 



\lxly)ly (Jx^z/^z ^j/(^j/^:e) ^z^z^x^) 

(iHl z ~\- lylx)^y ( iKly ~\~ Iz^xJ^z ^j/ ( itil z ^x^y) l z {ittly ~{~ Ix^z^) 

= 



[Li 2 ] = [i y ,i 2 ] = [i z ,P] = o 

From the commutators we learn that it is possible to measure the absolute 
value of the angular momentum and the component along one special di- 
rection at the same time. However, it is not possible to measure all three 
components at the same time. 

Central field problems: The potential V(r) is spherical symmetric and 
only depends on the distance from the origin. In this case we have 

[i x ,f] = [iy,f] = [Lf] = o 

[i x ,V(r)] = [iy,V(r)] = [i z ,V(r)] = 

[i x ,H} = [iy,H} = [i z ,H} = o 

[l 2 ,H]=0 



For central field problems H, I 2 and l z have a common set of eigenfunctions. 
Every stationary state is characterized by unique expectation values for H, I 2 
and l z . 



Angular momentum and Hamilton operators in spherical coordi- 
nates. 

Central field problems are best treated in spherical coordinates (R, 6, <fi). The 
operators in these coordinates are: 



(x,y,z) — > (R,9,(f>) . 
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We have to consider the rules of calculus, especially the chain rule. 
Results (without deduction): 



H = 



hd_ 

i d(f) 

-h 2 

2m 



1 d ( n d\ 1 d 2 1 
sm9— ) + 



sin 9 89 



d 2 2 d 



09 J ' sin 2 9 3(f) 2 

i 2 



+ 



2mr ; 



+ V(r) 



Therefore 



h 

I 2 > is a function of < 9, 
Hi I r,9,(f) 



Common set of eigenfunctions ty(r,9,4>). 



}' 



9 V independent part from eigenvalue problem for ^ I 2 



tf(r,0,<£) = R(r) 0(0) $(0) 



H 



H 



Eigenvalue problem for the operator l z 

% o<p 

Eigenfunctions for $(0) can be found easily: 



$(0) = Ce l 

h. 



= -imCe im<t> = mh$(<f>) 
i 

Initial conditions: 

$(0) = $(0 + 2tt) 



constants C, m 
eigenvalue mh 

periodicity 
m integer 
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Normalization: 

»2tt 



$*$ci0 = C' 







r-2-K 

I 



r-2-K 

= c 2 , 

Jo 



d(f) = 2nC 2 = 1 



C = (27T)- 1 / 2 

$(0) = (27i)- 1/2 e im ' t ' 



The z-component of the angular momentum can take the values mh (m in- 
teger) . 

Eigenvalue problem for the operator I 2 

Eigenfunctions Y{6, 0) and eigenvalues \h 2 



l 2 Y(9, 0) 



-h 2 

,2 



1 d 



sm6 09 V Sm 36 



d 



1 d 



2 1 



sin 6 d(j) 2 



Y{6A) 



\h 2 Y(9,(f) 



In the product function the eigenfunctions $(0) of l z are known. We get: 

d 2 



90 : 



-$(0) = -m 2 $(0) 



-h 2 



mnOdO \ oe 



+ 



d 2 



sin 9 d(j) 2 



e(^)$(0) = A/i 2 e(^)$(0) 



The operator is no longer dependent on and therefore we can divide by 
$(0). We get a differential equation for 0. This equation can be solved 
using similar methods as used for the harmonic oscillator. 

• Substitution x = cos 9 

• Power series ansatz 

• Recursion formula 

• Investigation of convergence, i.e. are there divergencies 

• Divergencies can only be avoided at special values of A where we get a 
finite power series. 

We get X — 1(1 + 1) with integer / > and / > \m\. The square of the 
absolute value of the angular momentum can only take values /(/ + l)h 2 with 
(I = 0,1,2,3,... and I > \m\). 
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Angular momentum eigenfunctions l^ m (^, (/>) 
Spherical harmonics Y l m (9,(f)) in product form: 

YT(M) = e(#)$(0) 

The 9 dependent part of 0(0) is proportional to the Legendre polynomials 
Pi for m = and the associated Legendre functions P ; ' m ' for m ^ 0. 
We get: 



X = cos 



ff"M = (W)' 



/7l m l 



Eigenvalue equation: 

zV(m) = i(i + i)H 2 Yr(9,d>) 

l z Y l m (9, ( p) = mhYr(9, ( p) 



(z = o, i, . . .) 

(|m| =0,1,2,...,/) 



(Z = 0,1,...) 

(-/ < m < Z) 



Spherical harmonics 

>T(M) = 



2/ + 1 (i- \m\y: 



An (l + \m\)\ 
Some special spherical harmonics: 

(if 



1/2 



pl m| (cos0)e^ 



^0° = I — I 



Y° 

1 2 



3 V /2 fl 

— cos 9 
An J 

3 \ 1/2 

— ) sinfle^ 
8n J 

5 \ 1/2 

' J (3cos 2 #-l) 



16tt 



/15\ 



1/2 



K^ 1 = ( — I sin^ cos 9e ±H 



Y 



±2 



\8n J 
15 \ 1/2 



32tt 



sm' 



p 
d 

d 

d 
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Linear combinations of complex spherical harmonics Y™ give the real spher- 
ical harmonics . 



il/2 i Y l 



The pair Y™, Y l m differs only in the factor e 
formulas: 



im0 p—imq. 



cos md> = —= \e~ 



sin rruj> 



iV2 



-im<j> g*" 1 ' 



we get using the definition of Y™ the real functions 



S/ m) (M) = N lm P{ mi (cos 9) 



M, 



cos m<; 
sinmd 



m > 
m < 

. Using the Euler 

m > 
m < 



m > 
m < 



General ansatz: 



H{r,0,<l>)*{r,9,<i>)=E*{r,0,<j>) 

9(r,6,<l>)=R(r)Yr(0,<l>) 

The angular part of the eigenf unctions is given by the spherical harmonics 
Y l m (9 : 4>) . With this ansatz and using the eigenvalue equation for r we get 
the radial Schrodinger equation for the hydrogen atom. 



'2fi 



d 2 2d 
dr 2 r dr 



](] + l)fr 2 Ze 2 
Rf r ) + HlZ ^ ig( r ) _ _i?( r ) = £i?(r) 
2/zr 2 r 



Solution of the radial Schrodinger equation 

• Using a new radial function u(r) = rR(r) results in the new equation 



h d 2 u(r) 
2/i dr 2 



+ 



1(1 + l)h 2 _ Ze 2 
2/ir 2 r 



u[r) = Eu(r) 



• Investigation of the asymptotic behavior for r — > oo and r — > leads 
to the following form of the solution 



u(r) = (Kr) l+1 e Kr w(Kr) 



k = 



-1\iE 



ft 2 



k is real for E < 0. 
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Inserting this definition into the radial Schrodinger equation leads to 
the Laguerre differential equation for w{kt). 



• Solutions of this differential equation can be found with a power series 
(recursion formula). 

• The power series is not diverging only for special values of E: Quanti- 
zation of energy. 

Ze 2 n , , i 

k + I + 1 = n 



2\E\ 

k, I, n integer; k > 0, / > 0, n > I 
Solution for E 



__zw 



n > 1 integer, n > I 
Energies 



Z 2 e 2 1 

E n = ~- , (n = l,2,3,...) 



2a n 2 

h 2 



a = — 

fxe z 



Bohr Radius = 0.529 A 



There is an infinite number of bound states with E n < 0. 
E = corresponds to proton and electron infinitely separated. 
E > non-bound states with additional kinetic energy. 
Ground state of hydrogen atom (Z — l,n = 1): 

e 2 

E x = w -13.6 eV 

2a 

Excited states of hydrogen atom (Z — 1, E n oc — ^): 

E 2 w -3.4eV, £ 3 ~ -1.5eV, ^ 4 ~ -0.85eV 



Each state ^ n ; m is characterized by three quantum numbers, however, the 
energy depends only on n. 
Energies are n 2 times degenerate. 
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n 


I 


m 


Notation 


States 


Degeneracy 


1 








Is 


1 


1 


2 








2s 


1 






1 


-1,0,1 


2p 


3 


4 


3 








3s 


1 






1 


-1,0,1 


3p 


3 






2 


-2,-1,0,1,2 


3d 


5 
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1.3.3 Radial function of hydrogen atom 



Rni(r) = - 



n-l-l)\ fIZ 



2n[(n + /)!] 3 \na 
Laguerre polynomials 



1/2 



\na I + V na 



dP 



Generalized Laguerre polynomials 



Explicite formulas 



Zr 

P = — 

na 



R 
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Rm 
R31 
-R32 



a J 



, 7 x 3/2 

i?2o = 2- 1 / 2 ^J (l-p)e-" 



9^ V a / 



9^ \ a ) 



3/2 



Is 
2s 
2p 
3s 
3p 
3d 
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Electron density distribution: 

Probability to find an electron within the volume dr: 

^* nlm y rdm dr = Kim(r, 9, 0)*„/ m (r, 9, <p)r 2 sm9drd9d<P 

Probability to find an electron within r and r + dr, independent of the angles 

(9,<f>): 

Pm{r)dr = R* nl R nl r 2 dr J J 4)^(0, 0) sin0d0# 

The integral is one (Normalization). For real R n i one gets for the radial 
probability density 

P nl {r) = R 2 nl r 2 

Probability to find an electron within the space angle dVt (within 9 and 9 + d9, 
(f> and (f) + d<f>), independent of the distance (r): 

W nl (0,<P)dn = Yr(0,<P)Yr(d,<P)dn J R n i(r)r 2 dr 

The integral is one (Normalization). With the well known formulas for the 
spherical harmonics Y™ we get: 

W nl {9,4>) = N 2 m [P! ml (cos9)] 2 = W lm {9) 

The angular distribution of the probability density is for the eigenfunctions 
of the hydrogen atom independent of <fi and therefore rotational symmetric 
w.r.t. the z axis. 




II = 2 
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Unit 



values in a.u. 



values in Si-units 



Charge elemental charge e = 1 e= 1.602189 • 1CT 19 C 

Mass electron mass m e — 1 m e =9.109534 • 10~ 31 kg 

Angular momentum Planck constant H H—l h= 1.054589 ■ 10~ 34 Js 

Length Bohr radius a = ^ = 1 a = 0.529177 • 10~ 10 m 

Energy 2\Ei(H)\ with fi = m e Hartree= ^ = 1 Hartree = 4.359814 • 10~ 18 J 



20 



Lecture 2 

Basic Mathematical Review 



2.1 Linear algebra 
2.1.1 Vector spaces 

A n- dimensional vector can be represented by its components w.r.t. a set of 
orthogonal and normalised basis vectors {e*j} 



a = e±ai + e 2 a 2 + h e n a n = ^ 



i=l 



For a given basis a vector is specified by its components. We represent a 
vector 

a 2 



— * 

as a column matrix. The scalar product of two vectors a and b is defined as 

n 

a ■ b = J^a*6j 



i=l 



This is due to e*j • = 5^ where 



1 if « = j 
otherwise 
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2.1.2 Operators 

An operator acts on a vector and transforms it into another vector 

Oa = b 

An operator is linear if 

O (xa + yPj = xOa + yOb 

A linear operator is completely determined if its effect on every possible 
vector is known. Because any vector can be expressed in terms of the basis, 
it is sufficient to know what O does to the basis vectors. Oe^ is a vector and 
can be expanded in the basis 



Oe, 



J2 

3=1 



1,2, 



. n 



The numbers Oji form a n x n matrix 



(O 



o = 



11 O12 

021 022 



o ln \ 

o 2n 



nl 



O 



n2 



O r 



O is the matrix representation of the operator O in the basis {e*j}. 

If A and B are the matrix representations of the operators A and £>, the 

matrix representation of the operator C = AB can be found as follows: 



C€j — ^ ^ CiCjj 



so that 



a 



i=i 

ABej 



A^e k B 

k=i 



i,k=l 



k=i 
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which is the definition of matrix multiplication and hence 



C = AB 



The order in which two operators or two matrices are multiplied is crucial. 
In general AB ^ BA and AB 7^ BA. This means two operators do not 
necessarily commute. The commutator of two operators is defined by 

[A, B] = AB - BA 

and can easily also be written in the matrix representation 

[A, B] = AB - BA 



2.1.3 Matrices 

A general rectangular n x m Matrix A has n rows and m columns 



Mi A 12 
A21 A22 



Al m \ 
A2m 



\A nl 



A 



n2 



A. r 



If n — m the matrix is square. Two matrices A and B can be multiplied 
C = AB if the number of columns of A is equal to the number of rows of B. 



AikB, 



kj 



k=l 



Matrix- vector multiplication is a special case of matrix-matrix multiplication. 

Some important definitions and properties of matrices. 

The adjoint of an N x M matrix A, denoted by A''", is a M x N matrix with 

elements 



fAt), 



A* 



i.e. we take the complex conjugate of each of the matrix elements of A and 
interchange rows and columns. If the elements of A are real, then the adjoint 
of A is called the transpose of A. 

The adjoint of the product of two matrices, is equal to the product of the 
adjoint matrices in reversed order. 



(AB)" 1 " = B f A f 
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Proof 



(AB)J, = (AB)* 



k k 

= E4X' = ( BtA %- 



The adjoint of a column matrix (vector) is a row matrix containing the 
complex conjugates of the elements. 

a = [a 1 , a 2 , • • • , a m ) 
If a and b are two column matrices with m elements, then 



i f b = «a£, ••• ,a* m ) 



\b m J 



i=i 



This is the scalar product of two vectors. 
Special properties of square matrices 

• A matrix A is diagonal if all its off-diagonal elements are zero 

A.. — A-fi- 

• The trace of A is the sum of its diagonal elements 

tr A = A a ■ 

i 

• The unit matrix is defined by = Sij and 

1A = Al = A 

• The inverse of A, denoted A -1 is a matrix such that 

A' 1 A = AA 1 = 1 

• A unitary matrix A is one whose inverse is its adjoint, a real unitary 
matrix is called orthogonal 



A" 1 = At 
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• A Hermitian matrix is self-adjoint, a real Hermitian matrix is called 
symmetric 

A f = A 

Further properties 

• trAB = trBA 

• (AB)- 1 = B X A 1 

• If U is unitary and B = U f AU, then A = UBU* 

• If A, B, C Hermitian matrices, and C = AB, then A and B commute. 

2.1.4 Determinants 

The determinant of an n x n matrix A is a number obtained from 
det(A) = | A 

where V{ is a permutation operator that permutes the column indices 1, 2, 3, . . . , n 
and the sum runs over all n\ permutations of the indices; pi is the number of 
transpositions required to restore a given permutation to natural order. 
Important properties of determinants 

• If each element in a row or in a column is zero the value of the deter- 
minant is zero. 

• If (A)ij = AuSij, then | A| = Mi = A n A 22 ■ ■ • A nn . 

• A single interchange of any two rows (or columns) of a determinant 
changes its sign. 

• |A| = (|At|)* 

• |AB| = |A||B| 

• If any two rows (or columns) of a determinant are equal, the value of 
the determinant is zero. 

. lA-^flAI)" 1 

• If AA t = 1, then |A|(|A|)* = 1. 

• If U f OU = D, and U f U = UU f = 1, then |0| = 



A ln 



Y J {-^TV l A ll A 22 ---A r , 
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2.1.5 Dirac notation 



Dirac introduced a notation for vector spaces that allows to express results 
in a concise and simple manner. We consider n basis vectors denoted by the 
symbol = 1,2, ... ,n, which are called ket vectors or simply kets. We 
assume this basis is complete so that any ket vector \a) can be written as 

n 

\°) = ^2 1*) ai ■ 

The column matrix a is the matrix representation of the abstract vector \a) 
in the basis We introduce an abstract bra vector (a\ whose matrix 

representation is a^. The scaler product between a bra (a\ and a ket \b) is 
defined as 



(a||6) = {a\b) = a f b = « a£, . . . , a* 
The square of the length of a vector 

n n 

(a\a) = J2 a *i a i = ^2 ' 0i 



b 2 
\b n ) 



i=l 



i=l 



i=l 



is always real and positive. In analogy to the ket basis we can introduce 
the bra basis {{i\} that is complete in the sense that any bra vector can be 
written as a linear combination of the bra basis vectors. 

(a\ = J2 a i (*l • 

i 

The scalar product now becomes 

(a\b) = J2 a i(i\j)bj ■ 



For this to be identical with the previous definition we must have that 
which is a statement of the orthonormality of the basis. 
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For a given ket \a) or bra (a\ the components of its matrix representation 
with respect to a basis {\i} or {{i\} can be determined by multiplication from 
the left or right with a general basis vector \j) or (j\ 

(j I a) = 5^0' M) °i = XI ^ ai = 

i i 

(« I i> = X < (i \j) = J2 a * 5i i = a ) 

i i 

Using these results we can write 

\a) =^\i) Oi = ^\i){i | a) 

i i 

and 

(a\ = J2 a i (*l = I 

i i 

which suggests we write 

j 

which is a statement of the completeness of the basis. 

An operator O is defined as an entity which when acting on a ket \a) converts 
it into a ket \b) 

0\a) = 1 6) . 

The operator is completely determined by its action on all basis vectors 

0\i) = Y,0 3t \j) = ^°h\j) 

3 3 

so that O is the matrix representation of the operator O in the basis {\i}}- 
The matrix representation is found by multiplying on the left by (k\ 

(k | O | i) = J> I j){Oh = EMO),* = (0) w 
which provides a useful expression for the matrix elements of O. 



2.1.6 Change of basis 

We specify two different basis sets, one, {\i)} where we use latin letters 
i, j, k, . . . for the vectors and another {\a)} where we use Greek letters a, f3, 7, . 
Thus we have 

(i\j) = $ij, ^2\i){i\ = l 
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and 

(a\P) = 5 a/3 , E \ a )( a \ = 1 • 

a 

Since the basis is complete we can express any ket in the basis {]«)} as 
a linear combination of kets in the basis and vice versa. That is, 

\a) = l\a) = J2 I a ) = E \ l ) U ™ = E NH U )*<* 

i i i 

where we have defined the elements of a transformation matrix U as 

(i\a) = U ia = (U)ia . 
Transforming in the opposite direction, we have 

i<> = i\t) = j2 1«><« i i) = E i")^ = E i«)( ut )- 

where we have used 

(a | i) = (i | a)* = W ia = (Ut) m . 

As a consequence of the orthonormality of the bases the transformation ma- 
trix U is unitary 

hi = (i I j) = E^ I a ^ a I 

a 

= Ef / -(u t W = (uu%. 

a 

which is in matrix notation 

1 = UU f . 

Starting from 8 a/ 3 — (a \ fl) one can show that 

1 = U f U 

and hence U is unitary. Thus we arrive at the important result that two 
orthonormal bases are related by a unitary matrix. The elements of the 
transformation matrix U are scalar products between the two bases. 
We now investigate the matrix representations of an operator for the two 
different bases. Suppose O is the matrix representation of O in the basis 
while fl is its matrix representation in the basis {]«)} 

Oij = (i\ O | j) ; n a p = (a\0\/3) . 
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We find the relation between O and ft from 



n a p = (a\0\P) = (a\lOl\P) 

= 5>|i>(i|0| J>O'I0> 

ij 

= ^(Ut) Qi (0)„(U)^ . 

ij 



Thus 



fl = u f ou 

or, multiplying on the left by U and on the right by 

o = umj f . 

These equations show that the matrices O and f2 are related by a unitary 
transformation. The importance lies in the fact that for any operator with a 
Hermitian matrix representation which is not diagonal in a given basis 
it is always possible to find a basis {|«)} in which the matrix representation 
of the operator is diagonal, i.e. 

= Ua8aj3 ■ 



2.1.7 Eigenvalue problem 

If the result of the operation O | a) is simply a constant times \a), i.e. 

O | a) — uj a \a) 

then we say that \a) is an eigenvector of the operator O with an eigenvalue u a . 
Without loss of generality we can choose the eigenvectors to be normalised 

(a | a) = 1 

• The eigenvalues of a Hermitian operator are real. 
This follows from 

(a | O | a) = {a \ \ a) = {a \ O | a)* 

Multiplying the eigenvalue relation by (a\ and using the above relations 
we have 

which is the required result. 
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The eigenvectors of a Hermitian operator are orthogonal. 
We consider 

and its adjoint 



<J3\0* 



Since O is Hermitian and ujp is real, we have 

((3 | O = {fi\ufi 

Multiplying O \ a) — co a \a) by (f3\ and ((3 \ O = (f3\uop by |a) and 
subtracting the results we obtain 

(u/3-u a )(P I a) = 

so that (/3 | a) = if w Q ^ up. Thus orthogonality follows imme- 
diately if the two eigenvalues are not the same (nondegenerate). It 
can be shown, that degenerate eigenvectors can always be chosen to be 
orthogonal (we omit this here). 

The eigenvalue problem can be posed as follows. Given the matrix represen- 
tation, O, of a Hermitian operator O in the orthonormal basis we wish 
to find the orthonormal basis {\a)} in which the matrix representation, fl of 
O is diagonal. In other words we wish to diagonalise the matrix O. The two 
representations of the operator O are related by an unitary transformation 

fl = u f ou 

Thus the problem of diagonalising the Hermitian matrix O is equivalent to 
the problem of finding the unitary matrix u that converts O into a diagonal 
matrix 

u \ 

u 2 

U f OU = cu = 

o ••. 

The problem can be reformulated as follows. Given the n x n Hermitian 
matrix O, find all n distinct column vectors c (the eigenvectors) and corre- 
sponding n numbers u (the eigenvalues) such that 

Oc = LUC . 

This equation can be rewritten as 

(O-cj1)c = . 
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Nontrivial solutions (c ^ 0) only exist when 

|0 = . 

This determinant is a polynomial of degree n in the unknown uj. It has n 
roots uj a which in this case will be the eigenvalues of the matrix O. Once 
we have found the eigenvalues, we can find the corresponding eigenvectors 
by substituting each oo a into the above equation and solving the resulting 
equation for c a . In this way c a can be found within a multiplicative constant, 
which is finally determined by requiring the eigenvectors to be normalised 

J2«y< = 1 • 

i 

In this way we find n solutions to 

Oc a = uo a c a a — 1, 2, . . . , n 

Since O is Hermitian, the eigenvalues are real and the eigenvectors are or- 
thogonal. We construct a matrix U defined as U ia = cf , i.e. 







c\ ■ 


■ cf\ 


u = 


4 


c\ ■ 












{4 


cl ■ 





Thus the ath colum of U is just the column matrix c a . Then from the 
orthogonality relation we find 

E^^ = E( ut )«( u )^ = ^ 

i i 

which in matrix notation is 

U f U = 1 

Multiplying OU = Uoo by U f we have 

U f OU = u 

which shows that U is the unitary transformation between the original and 
new basis for the diagonal representation of O. 
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2.1.8 Functions of matrices 



Given a Hermitian matrix A, we can define a function of A, i.e. /(A), in 
much the same way we define a function f(x) of a simple variable x. For 
example, the square root of a matrix A, which we denote by A 1 / 2 , is simply 
that matrix which when multiplied by itself gives A 

A l/2 A l/2 = A 

The sine or the exponential of a matrix are defined by the Taylor series of 
the function 

exp(A) = 1 + -A + -A 2 + -A 3 + • • • 



or in general 



/(A) = ]T Cfc A* 



k=0 



We are still faced with the problem of calculating A 1 / 2 or exp(A). If A is a 
diagonal matrix 

(A)jj = Oj<5jj 

everything is simple, since 



(a\ 



\ 



a\ 

•. 



so that 



/E 



k=0 



k c k a\ 



V 



Efc Ck(4 o 
o 



f(a 2 ) 




\ 



\ 



f(a n )J 



Similarly, the square root of a diagonal matrix is 

A x / 2 = 



af 







V 



1/2 , 
an ) 



32 



If A is Hermitian we can always find a unitary transformation that diago- 
nalises it 

U f AU = a 

The reverse transformation is 

A = UaU f 



Now notice that 
or in general 
so that 



A 2 = UaUtUaUt = Ua 2 Ut 
A k = Ua^ 



/(A) = c k A k = U PT c fe a fc Ut = U /(a) U* 



= U 



k=0 \ k 

tfM \ 

f(a 2 ) 


V fMJ 



Thus to calculate any function of a Hermitian matrix A, we first diagonalise 
A to obtain a, the diagonal matrix containing all the eigenvalues of A. We 
then calculate / (a), which is easy because a is diagonal. Finally we transform 
/(a) back using the matrix U formed by the eigenvectors of A. For example, 
we can find the square root of a matrix A as 

A l/2 = Ua l/2ut 

since 

A l/2 A l/2 = Ua V2 U t Ua l/2 U t = Ua l/2 a l/2ut = UaU t = A 

2.2 Perturbation theory 



For an un-perturbed model system with the Hamiltonian the exact wave 

rgies En^ are knowi 



functions and exact energies En^ are known. 
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For the system of interest, a perturbed system with the Hamiltonian H, the 
wave functions \l/ n and energies E n cannot be calculated exactly. 



Perturbation theory describes how the exact solutions of the model system 
change under the influence of a perturbation (perturbation parameter 

A, perturbation operator H^>). In a systematic procedure corrections to the 
energy and the wave functions in first, second , third, . . . order are defined: 

eV,eV,E®,... and *« *W *W ... 

The perturbation XH^ has to be small and we have 

H = + 

E n — En* 1 _^ 

2.2.1 Derivation 

Ansatz: 

E = E<® + A_E (1) + \ 2 E^ + ■■■ 
* = #(°) + A^ (1) + A 2 ^ (2) + ■■■ 
= (*(°) | withp= 1,2,... 

Formal power series in the perturbation parameter A. Orders of magnitude: 

|E(°)|»|A J E( 1 )|»|A 2 E( 2 )|»... 

Insert in the Schrodinger equation: 

(if (°) + Aif + A*« + A 2 ^ 2 ) + •••) = 

( S (0) + XE (1) + A 2 E (2) + . . . ) ^(0) + A ^(l) + A 2^(2) + . . . ) 

Combine w.r.t. the orders of powers of A: 
"#(<%(o) _ ^(0)^(0)1 + 

A [fl-(o)^(i) + ^-(i)^(o) _ £(0)^(1) _ ^(1)^(0) 



+ A 2 



#(0)^(2) + #(1)^(1) _ #(0)^(2) 

_ £(2)^(0)] + ... = 
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For an arbitrary perturbation 
#((%(°) 

+ ^(i)^(o) 

jf (0)^(2) + ^(1)^(1) 

+ ^(l)^(P-l) = £(<%(p) + ^(l)^(p-l) _| |_ E (p)yj/(0) 

The solution for the un-perturbed problem is known. For the calculation of 
perturbative corrections the equations above have to be solved. 
Solutions: Equations are multiplied from left with and then integrated 

over all coordinates. 

2.2.2 Energy corrections 

First order: 

The right hand side can be simplified using that the unperturbed wave func- 
tion \I/( ) is normalised and orthogonal to the corrections vj/fr). The first term 
on the left hand side is zero because is Hermitian and the orthogonality 
condition. 

| if(°)tf(i)) = (jy(o)^(o) | = ^(0)^(0) | = o 
Therefore we get 

E W = ^(0) | ^(i)^(o)^ 

The energy correction in first order is equal to the expectation value of the 
perturbation operator w.r.t. the un-perturbed wave functions. 
Higher orders: 

From using the same derivation we get: 

E {2) = ^(0) I #(1)^(1)) 
E (p) = ^(0) I #(l)^(p-l)) 

The wavefunction in (p — l)-th order determines the energy in p-th order 
within perturbation theory. 

Note: A more thorough analysis shows that the wave function in p-th order 
is sufficient to calculate the energy up to (2p + l)-th order. 



parameter A each order has to be zero. 

=£(<%(2) + ^(1)^(1) + ^(2)^(0) 
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2.2.3 Corrections and E { n ] 

y&n^ is expanded in a complete set of solutions to the un-perturbed problem 

^(0)_ 

With this orthogonality of \E , i°' ) w.r.t. ^i 1 ^ is guaranteed. 

« 0) K 1) > = E a i« 0) i*f> = o 

From the general equation in first order we get: 

^(1)^(0) _ ^(1)^(0) = £(0)^(1) _ #(0)^(1) 
n n n n n n 

= (4°» -i/ (o ')i>*r 

= E«,K ) -^ 0, )*f 

Now we use the usual procedure: Multiplication from left with (i ^ n) 
and integration. 

(*(«» | ^(1)^(0)) = at ^(0) _ = H w 

11 

definition we get: 



This defines the matrix element of the perturbation operator. With this 



(i) 



P (0) p (0) 



Energy correction in second order: 



o-(i) 

$W = V ^ ^ (0) 



rr(l) 



p(0) _ p (0) 



I in I 



P (0) _ p (0) 
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2.2.4 Summary 

The equations of perturbation theory contain the energies E [ n ] and wave 
functions of the un-perturbed problem and the matrix elements of 
the perturbation operator: 

= <*<°> | HW*W) 
Perturbation series of the energy: 



E = E (0) + H {1) + V 1 in 1 + 

Cjn ^ + fl ™ + ^(o) _ ^(o) + 



Perturbation series of the wave function: 



*n *n I" ^ (o) _ (0) *t + 

Higher order corrections are not used very often. 



2.3 The variation method 

The expectation value of the energy E^ calculated with an arbitrary (valid) 
wave function \1> is an upper bound for the exact energy E of the ground 
state of the system. 

Ey = | H^) > E * normalised 

Proof: Expand ^ in the complete set of the exact eigenfunctions \I/ n of the 
system. 

H^b n E n ty n 

oo 

# = a n^n 
n=0 

Normalisation: 

oo oo 

n=0 m=0 

oo 

(*i*) = Ew 2 = 1 

n=0 
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Expectation value of the energy: 



n=0 m=0 

oo oo 

= J2J2 a n amEm ^ n I ^ 

n=0 m=0 

oo 

— ^2 i a «i 2 -^ n 



n=0 

In the last equation we replace all eigenvalues E n with the lowest eigenvalue 
E (ground state), with E < E 1 < E 2 Then we have 

oo 

Eq, > Eq | a n | 2 = Eq 

n=0 

2.3.1 The variation principle 

Within all valid wave functions \I/ for the ground state we are looking for the 
function with the lowest expectation value for the energy Eq,. 

Eq, — > Minimum 

Practical procedure: Choose an Ansatz for ^ with free parameters ctj and 
determine the parameters using the variational principle. This gives the best 
possible solution within this Ansatz. 

This is equivalent to a extremal value problem. The variation SEq, has to be 
zero for any arbitrary variation of the parameters ctj. 

r 7-i / \ dEq, dEq, 

dEq,(a 1 , a 2 , ■ ■ ■) — t; — ««i + t; — da 2 H = 

ooii oa 2 

Because the variations dai are arbitrary, we have: 

= for all % 

OOL; 



2.3.2 Linear variations 

Variational function \1/ is chosen as a linear combination of n fixed basis 
functions $j: 

n 

i=i 
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The coefficients Cj are varied until the minimal value of the energy expectation 
value is found. For not normalised functions \& we have: 



E 



| *) B 

fc j 



Notation: 

iJ fci = ($ fc | =#* fc (Matrix element of #) 

S'fcj = | $j) =S'* fc (Overlap integral) 

Minimisation of E<$ w.r.t. the variational parameters q: 
dEv 



do 



for i — 1, 2, . . . , n 



dEy 1 /'Si A dB \ 1 P A A<9B> ] 



<9cj B 2 \ dci dci J B \dci B dci 
Because B ^ (Normalisation integral) and with Eq = A/B we get: 

dA 8B _ 

o — (J 

OCi OCi 

Calculation of the derivatives: 

k j 

i 

B = ^2J2 c *k c i Sk i 

k j 

Inserting the derivatives we get the generalised eigenvalue equation: 
(ifjj — E^Sij) Cj = for i = 1, 2, . . . , n 
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Known: Integrals and 

Unknown: Coefficients Cj and corresponding energies 

Linear algebra: Non-trivial solution of a homogeneous system of linear equa- 
tions, for the case of vanishing determinant: 

\Hij — Eq,Sij\ = 

The solution of this equation gives n energies E k , where the lowest value is 
the variational energy of the ground state. For each energy E k we get after 
inserting in the system of linear equations the coefficients Cj k and therefore 
the corresponding wave function ^ k . 
Additional index to distinguish the individual solutions: 

n 

^ HijCij = ^2 SijC jk E k for i,k = 1,2,..., n 

j j 

Using matrices: 

H = {Hij}, S = {Sij}, C = {Cij}, E = {Ei5ij} 
HC = SCE 

Using the nxn -matrices H, C, S and E we can write the equation in matrix 
form. H and S are real for real basis functions $j and symmetric (hermitian 
otherwise), E is a real diagonal matrix. 
Important special case: Orthonormal basis 

HC = CE 

H is the matrix representation of the Hamilton operator in the orthonormal 
basis. Through diagonalisation of H we get the energy eigenvalues (diagonal 
elements of E) and eigenvectors (columns of C). We have 

E = C~ l HC = C T HC 

n 

C T C = I or^|c ifc | 2 = l for k = 1,2,..., n 

i=i 

The inverse and transpose of the matrix C are identical. Eigenvectors are 
orthonormal. 
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Lecture 3 



Molecular Hamiltonian 



3.1 Time-independent non-relativistic Schro- 
dinger equation 

For a system with N electrons and K nuclei the time-independent Schrodinger 
equation is 

H(r, R)V n (r, s, R, a) = E n ^ n (f, s, R, 3) . 

Notation: 

r*j = (xi, i/i, Zi) Coordinates for electron i 

Si Spin - coordinate for electron i 

Ra = (%a, HA, za) Coordinates for nuclei A 

a a Spin - coordinate for nuclei A 

r — (r*i, . . . , fjv) Coordinates for all electrons 

s— (si, . . . , sjv) Spin - coordinates for all electrons 

— * — * — * 

R = (Ra, • • • , Rk) Coordinates for all nuclei 
= ipAi ■ ■ ■ i &k) Spin - coordinates for all nuclei 

Conventions: 

Electrons : Number iV ; Indices . . .; Mass m 
Nuclei: Number K ; Indices A,B, . . .; Mass Ma 

Non-relativistic Hamilton operator 

H(f, R) = f N (R) + f cl (f) + V Ne (f, R) + V ee (f) + Vn N (^) 
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Kinetic energy: 

Nl } " " V 2M a \dx\ + 8y\ + dz\) " ' ^ 2M A V 

l ^ ~ 2m^\dx 2 l + dyf + dzf ) 2m ^ 1 

Potential energy: 

A i 



e 2 



■ 3" ^ 



e 2 



Ra r 

A<B AU 

Where Za is the charge of nuclei A and = |rj — Tj\. 

H includes the non-relativistic kinetic energy and the potential energy of the 

Coulomb interaction. 



Not included are: 

• Relativistic correction to the kinetic energy 

• Interaction of magnetic moments 
(orbit /orbit, Spin/orbit, Spin/Spin) 

• Interaction with external electric and magnetic fields 



3.2 Born— Oppenheimer Approximation 
3.2.1 Part I 

Goal: Separation of nuclear motion and electron motion. 

Idea: Electrons are much lighter than nuclei and move therefore much faster. 
They adjust to each new position of the nuclei almost instantaneously. 

Ansatz: The electronic wave function ^ e \(r; R) is calculated for fixed posi- 
tion (R) of the nuclei. It directly depends on r and parametrically on 
R. 
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Electronic Schrodinger equation 

H cl (f; R)* A (r; R) = E cl (R)^ el (r; R) 
with the electronic Hamilton operator 

H cl (f; R) = f c i(r) + V Ne (?; R) + V ec {7) . 
For the derivation we have used 

• T N (_R) = non- moving nuclei 

• Vnn(-R) does not act on r (multiplicative operator) and can be neglected 
for the solution 

Interpretation of the equations 

• The calculation of ^ci takes into account the electronic interaction of 
nuclei and electrons, however not their dynamic coupling. There is no 
exchange of kinetic energy. 

• Macroscopic analogy: No heat exchange, adiabaticity. 

• The solution of the electronic Schrodinger equation is difficult but pos- 
sible (Quantum chemistry). 

• The eigenfunctions ^(r; R) are a complete ortho normal set of func- 
tions. 

3.2.2 Part II 

The total wave function \l/(r, R) is expanded in the electronic wave functions 
*3(r; R) 

V(r,R) = Y,Xn(R)n(r;R) . 

n 

Insertion in the Schrodinger equation 

H(f, R) J2 Xn(R)n(r; R) = EJ2 Xn(R)n(r; R) 



Multiplication from left with \E'™* (r; R) and integration over all electronic 
coordinates r 

f N (R) + E™(R) + V m (R) - E] Xm(R) = ~ C mn (R) X n(R) , 
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and 



Note for the derivation: The coupling terms C mn (R) are due to the action of 
f N {R)on^(f;R). 

Comparison of relevant matrix elements shows that C mn (R) is normally van- 

ishingly small compared to f N (R) + E™(R) + Vnn(-R)] • 

Possible approximations 
Cmn = for m 7^ n : Adiabatic approximations 
C mn = for all to, n : Born-Oppenheimer approximation 

Both approximations lead to non-coupled systems of equations for the nuclear 

functions Xn(R)- 



3.2.3 Part III 

After the decoupling of the nuclear Schrodinger equation we get for each 
electronic state an effective Schrodinger equation for the nuclear motion: 

H N (R) X m(R) = EXm(R) , 

with 

H N (R) = f N (R) + V e s(R) , 
where we have been using the adiabatic approximation 

V cS (R) = E™(R) + Vn N (£) + C mm {R) 

and the Born-Oppenheimer approximation, respectively. 

V eS (R)=E™(R) + V NN (R) 

The diagonal correction term C mm (R) is small and difficult to calculate. Most 
calculations are done using the Born-Oppenheimer potential energy surface: 

U(R) = E™(R) + V N x(R) 

The electronic energy E™(R) and the interaction energy of the nuclei are 
independent of the nuclear masses, therefore we see that 

Isotopes have the same potential energy surface U(R). 

The Born-Oppenheimer approximation is often valid to high accuracy. Ex- 
ceptions occur if different electronic states are significantly coupled through 
nuclear motion. Special care is necessary for 
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• close lying potential energy surfaces (crossings and avoided crossings 
of electronically excited states) 

• rapid movements of nuclei (impact of fast particles). 

3.2.4 Order of magnitude of effects for H2 

a) Relativistic effects in H 2 (State X X E+) at the equilibrium geometry R e 

Change of D : -0.48 cm -1 

£>i : -0.34 cm" 1 

E(u = 1) - E(v = 0): -0.14 cm" 1 

b) Non-adiabatic effects in H 2 (State X l T^) at the equilibrium geometry 

Re 

Approximately calculated change in energy 

H 2 : 0.65 cm" 1 , ED : 0.42 cm" 1 , D 2 : 0.23 cm" 1 

All effects are certainly smaller than 1 cm -1 , 

Dq = 36118.3 ± 0.5 cm -1 (experiment) 

D = 36118.0 cm -1 (theory; adiabatic + relativistic) 

c) Born-Oppenheimer effects, more precisely: energy change due to the di- 

agonal adiabatic correction C nn 
H 2 X X E+ 120 cm" 1 
a 3 E 9 100 cm" 1 
5 3 E+ 200 cm" 1 
E, E 1 ^ 430 cm -1 (avoided crossing) 

This could be important for the exact calculation of electron spectra 
(80 cm" 1 ss O.OleV). 

Not important for vibrational spectroscopy, Born-Oppenheimer effects 
on calculated transition energies (Au = 1) are smaller than 1 cm" 1 for 
H 2 (State A^E+) for v = 1 - 10. 

3.3 Stationary points 

Atoms, 3A degrees of freedom (coordinates) 
For non-linear molecules: 

• 3 degrees of freedom for translation 

• 3 degrees of freedom for rotation 
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• 3iV — 6 degrees of freedom for internal coordinates 
Coordinates (e.g. Cartesian coordinates): 

R — {Rii R2, ■ ■ ■ , Rzn) = {Ri} 



Definition 



9* 

fij 



dU(R) 

dR l 
d 2 U(R) 
ORidRj 



components of gradient 
harmonic force constants 



Taylor expansion within the vicinity of a reference point with coordinates 
Ro = {Ri,o}'- 



U = U(0) + ^2 9i x i + g fn x i x i + 

with 

or in matrix form 



* 3 



x i — Ri ~ Ri,0 



U = U(0) + g T x + -x T Fx + 



with 









( 9 A 




x 2 




92 


X = 




9 = 






\X3N J 




\93nJ 



F = 



37V, 1 



Fl,3N 



Definition: 

Points with the coordinates R e = {Ri, e } are called stationary if all compo- 
nents of the gradient are zero. 

gi (R e )=0 fori = 1,2,..., 37V 

At a stationary point we have the Taylor expansion 

» 3 
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with 

Xi Rj, Ri,e • 

Characterisation of stationary points: 

The force constant matrix F(R e ) has 3N real eigenvalues, whereof 

• 3 for translation (equal zero) 

• 3 for rotation (equal zero) 

• 3iV — 6 for internal motion (different from zero) 

The number of negative eigenvalues is used to characterise the stationary 
point 

• for a minimum 

• 1 for a transition state 

• 2, 3, ... for a saddle point of higher order 

3.4 Internal coordinates 

Changes in bond length and in the angles between chemical bonds provide 
the most significant and physically meaningful set of coordinates for the 
description of the potential energy. These coordinates are called "internal 
coordinates" since they describe (being unaffected by translation and rota- 
tions of the molecule as a whole) just the internal motions of the molecule, 
i.e. the molecular vibrations. 

The types of internal coordinates which are generally used are the following: 

• Bond stretching coordinate : variation of the length of a chemical bond. 




• In-plane bending coordinate : variation of the angle between two chem- 
ical bonds having one atom in common. 
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• Out-of-plane bending coordinate : variation of the angle between the 
plane defined by two bonds with one atom in common and a third bond 



connected to the common atom. 




• Torsion coordinate : variation in the dihedral angle between the planes 
determined by three consecutive bonds connecting four atoms. 
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A 




Since the variations in bond length and angles in a molecule can be written 
in terms of the variation in the Cartesian displacement coordinates of the 
atoms, the problem we have to face is that of writing the explicit form of 
the transformation from Cartesian coordinates to the new set of internal 
coordinates. 

Let us call Sk a generic internal coordinate. The most general relation be- 
tween s k and the Cartesian displacement coordinates can be written in the 
form 

1 . 

Sk = B\ Xi + — BijXiXj + higher terms 

i ij 

where the coefficients Bf,Bfj, etc., are determined by the molecular geome- 
try. A drastic simplification can be achieved if we restrict our treatment to 
the case of infinitesimal amplitudes of vibration where we can drop all terms 
not linear in x. 

If we call s and x the vectors whose components are the internal and the 
Cartesian coordinates, the linear equation in matrix notation is 

s = Bx 

where B is a matrix whose elements are the coefficients B\. 
The matrix B is not, in general, a square matrix. Let us consider, for example 
a bent triatomic molecule and choosing as internal coordinates the two bond 
length and the bond angle. There are also 3n — 6 = 3 vibrational degrees 
of freedom and with the 9 Cartesian coordinates the B-matrix will have 3 
rows and 9 columns. The translations and rotations are not included in the 
set of internal coordinates which describe, by definition a molecular motion 
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in which the relative positions of the atoms are changed. Since B is not a 
square matrix, it cannot be inverted. 

We can always include in the s-vector six additional coordinates to describe 
the three translations and the three rotations. These six coordinates are 
most conveniently chosen in the form 




where M = Yli m i * s the total mass of the molecule and I x , I y , and I z are 
the moments of inertia. and are the coordinates of the equilibrium 

position of the 7th atom relative to the centre of mass and 

qix = \frn~ikxi q Iy = ^/raJA?// q Iz = y/mJAz! 

In this way B is a square and non-singular matrix (i.e. |B| 7^ 0) and therefore 
can be inverted. We can thus define the inverse transformation to get from 
internal to Cartesian coordinates 

x = B x s 

where B _1 is the inverse matrix of B including translation and rotation. 
3.4.1 Z-matrix 

A special choice of internal coordinates often used as input to quantum chem- 
ical programs is the Z-matrix. The Z-matrix of a molecule is build from a 
list of the atoms. The coordinates of each new atom (A) are characterised 
by 

1. The distance to another, already defined atom B 

2. The angle between the bond to atom B and the bond of atom B to an 
already defined atom C. 

3. The dihedral angle between A, B, C, and an already defined atom D. 
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The first three atoms are special cases with 

• The first atom is at the origin of the Cartesian coordinate system. 

• The second atom is defined by the distance to atom 1 and is oriented 
along a Cartesian coordinate axis. 

• The third atom is defined by the distance to atom 1 or 2 and the angle 
(3-2-1) or the angle (3-1-2) 

The Z-matrix of a molecule is not unique. There are many different possibil- 
ities to order the atoms and to define distances, angles and dihedrals. 

Example: Z-matrix for staggered conformation of ethane. 




1 


C 




















2 


C 


1 


54 


1 














3 


H 


1 


01 


1 


109 


5 


2 








4 


H 


1 


01 


2 


109 


5 


1 


180 





3 


5 


H 


1 


01 


1 


109 


5 


2 


60 





4 


6 


H 


1 


01 


2 


109 


5 


1 


-60 





5 


7 


H 


1 


01 


1 


109 


5 


2 


180 





6 


8 


H 


1 


01 


2 


109 


5 


1 


60 





7 



In the first line of the Z-matrix atom 1, a carbon atom, is defined. Atom 
number 2 is also a carbon that is a distance of 1.54 A from atom 1 (columns 
3 and 4). Atom 3 is a hydrogen atom that is bonded to atom 1 with a bond 
length of 1.01 A. The angle formed by atoms 2-1-3 is 109.5°, information that 
is specified in columns 5 and 6. The fourth atom is a hydrogen, a distance 
of 1.01 A from atom 2, the angle 4-2-1 is 109.5°, and the torsion angle for 
atoms 4-2-1-3 is 180°. Thus for all except the first three atoms, each atom has 
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three internal coordinates: the distance of the atom from one of the atoms 
previously defined, the angle formed by the atom and two of the previous 
atoms, and the torsion angle defined by the atom and three of the previous 
atoms. 

The Z-matrix can often be written in many different ways as there are many 
combinations of internal coordinates. As a rule internal coordinates should 
be chosen to have minimal coupling. However, it should be stressed that the 
final results of a quantum mechanical calculations are not dependent on the 
actual coordinate system chosen. 
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Lecture 4 

Two-Electron Systems and Spin 



4.1 Helium atom: Ground state 
4.1.1 Basics 

Electronic Hamilton operator within Born-Oppenheimer approximation: 

H cl (f; R) = f c (r) + V Ne (f; R) + V ee (r) 

Coordinate system: nuclei at origin, spherical coordinates for electron 1 
(r 1? 61, 0i) and electron 2 (r 2 , 6 2 , 2 ). 
Electronic Hamilton operator for He (Z=2): 

f, ,- -> fi' 2 m2 „2n Ze 2 Ze 2 e 2 
^•'' 2 » = ^< V ' + V2 »- — 'IT + ^ 

Corresponding Schrodinger equation: 

^el(n,f 2 )* e l(rL,f 2 ) = E e i^ c i(fi,f 2 ) 

For this differential equation with 6 variables there is no exact solution. We 
are looking for approximations, first for the ground state. 

Experiment: 

The energy E of the ground state is equal to the sum of ionisation potentials 
for He (24.59 eV) and He+ (54.42 eV). 

E (exp) = -79.01 eV 

This is our reference to assess different approximations. 
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4.1.2 Orbital model 



Ansatz: 

*ei(r!,r 2 ) =*i(fi)* 2 (r 2 ) 

One electron functions are called orbitals. In the orbital model the N-electron 
function is approximated by a product of N orbitals. 
Probability interpretation: 

^ cl dr = [^ 1 (f 1 )*^ 1 (n)dr 1 ] ^ 2 (f 2 )^ 2 (f 2 )rfr 2 ] 

w 12 = w x w 2 

The probability W12 to find at the same time electron 1 at f\ and electron 2 

at r 2 is equal to the product of the single probabilities W\ and Wi. 

In the orbital model the probability distribution of an electron is independent 

from the other electrons ('independent particle model'). 

Note: 

More accurate methods that go beyond the orbital model have to include the 
correlation of electrons. 



4.1.3 Formal separation 

The Hamilton operator consists of one- and two-electron parts 

H el = + V ee = H {0) + — 

# (0) = 5>(rl) 
i=i 

i h2 ^2 Ze 2 

k * = -^-- 

The one-electron operator hi is the same as the Hamilton operator of a 
hydrogen-like atom with Z = 2 for He. Eigenvalues and eigenfunctions \I/ * 
are known. 

h^kifi) = e fc * fc (fi) 
Z 2 e 2 1 

e k = - with m = 1,2,3, .. . 

2a n% 

ty k = Tr- 1/2 (Z/a) 3/2 e-^ for n k = 1 

Crude Approximation: With V cc = the Hamilton operator separates like 
H^. The two-electron wave function can in this case be written as a product 
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of orbitals (the solutions for the hydrogen atom with Z — 2), and we get for 
the energy 

£ (0) =ei(m) + e 2 (n 2 ) 
Ground state of the He atom with n\ — n 2 — 1: 

Z 2 e 2 

= = -108.84 eV 

(X 

7rcr 

Because of the assumption V^e = we get a huge discrepancy to the experi- 
mental ground state energy (—79.01 eV). 



4.1.4 Perturbation theory 

We take the electron-electron interaction V^, as the perturbation operator 
This is not really a valid choice as V^e is not small compared to the other 
terms in the Hamilton operator. Nevertheless we calculate the correction to 
the energy to first order: 



fl2 



Inserting the analytic solutions for leads to a difficult six-dimensional 
integral for the variables (n, 9i, fa, r 2 , 2 , fa), that can be solved analytically. 



5Ze 



2 



8 a 



For He with Z = 2 we get: 



= 34.01 eV 

+ = -74.83 eV 

With the perturbation correction to first order the energy is off by 5.3% from 
the experimental value (-79.01 eV). 

Corrections to higher order are possible, but much more complicated than 
simple variational calculations. 



4.1.5 Variational calculations 

Idea: In the solutions used up to know, both electrons were seeing the field 
of the nuclei with Z = 2. This neglects the partial screening of the 
nuclear charge due to the other electron. 



55 



Ansatz: We replace the nuclear charge Z = 2 in the wave function for the 
ground state by an effective nuclear charge rj < 2. 



7TGT 



Variational calculation for the determination of the optimal value of rj: 

dE^j d . . - 

— = -(*|// cl |*) = o 

The calculation of the expectation value of H c \ w.r.t. to \I/ affords the solution 

of several integrals (difficult but analytic solutions exist). 

Solution: 



77 — 2t]Z + -77 
8 



e 
a 



2^ - 2Z + - 



e 

— = 

a 



^o P t - Z - — 



Eopt = -V, 



opt" 



5 \V 



16 a 



For He with Z = 2 we get: 



77 op t(#e) = 1.6875 
E opt (He) = -77.49 eV 



The nuclear charge is screened by 0.3e. The error in the ground state energy 
is reduced from 5.3% with rj = 2 and a perturbative calculation to 1.9% with 
V = Vopt- 



4.1.6 Overview of ground state energies 

Goal: Assess accuracy of approximative solutions by comparison to experi- 
ment. 
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£j (a. u.J 


error 


Experiment 


-Z.yUo t o4:{6 J 




Orbital model 






Zeroth order (without V cc ), E^ 


-4.00 


-37.8 


Perturbation theory 1. order, + E^ 


-2.75 


5.3 


Simple variational calculation, E opt 


-2.8477 


1.9 


Best variational calculation 


-2.8617 


1.4 


Including electron correlation 






Simple variational calculation, 40 Parameters 


-2.90276 


0.04 


Best variational calculation, 1078 Parameters 


-2.9037462 


0.001 


With additional relativistic corrections 


-2.9037848 





4.2 Electron spin 

Experimental fact: 

Electrons have a proper angular momentum (Spin), the components of which 
can have only two distinct values in any arbitrary direction: 




Quantum mechanical description with spin operators: 

s 2 = §1 + Sy + s 2 = s+s_ - s z + s 2 z = s-s + + s z + s 2 

S_|_ S X -\~ ISy 

S_ S X ISy 

The linear combinations s + and s_ are useful for derivations. Spin operators 
have the same commutation relations as the angular momentum operators: 

[s 2 , s x ] = [s 2 , s y ] = [s 2 , s z ] = 
Eigenfunctions: 

For a one electron system there are two orthonormal eigenfunctions |a) and 
\(3). Application of operators: 
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\a) 




s 2 


fl«> 


m 


Sz 


\\a) 


-m 




m 


\\a) 


Sy 


m 


-\\a) 







\a) 


S- 








General eigenvalue equation for spin: 

s z \a) = m s h\a) 
§ 2 \a) = s(s + l)h 2 \a) 

Spin eigenfunction \a), quantum numbers s and m s for the total spin and its 
z-component with — s < m s < s. 
For a single electron we have: 

1 1 

m s — s = - 

2 2 

Spin orbitals: For the complete description of an electron we need three 
spatial variables and one spin variable. The spin orbitals \I/ (f) \a) are products 
of a spatial orbital *ff(r) and a spin eigenfunction \a), \a) oder \/3). 

Two electron systems: 

Operators and eigenvalue equation: 

S x = $xl + S X 2', S y = S y i + S y 2 ', S z = S z \ + S 2 2 

S 2 = S 2 + S 2 + S 2 Z 

Because we have the general commutator relations for angular moment op- 
erators, the eigenvalue equations have the form: 

S>i(j 2 ) = M s %i<7 2 ) (S<M S <S) 
S 2 \a 1 a 2 ) = S(S + l)H 2 \a 1 a 2 ) 
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Spin eigenfunction |o"i<7 2 ), quantum numbers S and M s for total spin and its 
z-component. 

Relation to the Hamilton operator: 

We use a non-relativistic Hamilton operator without any spin components, 
therefore Hamilton operator and spin operators commute. 

[H,S 2 ] = [H,S z ]=0 

Conclusions: 

• There are exact expectation values for S 2 and S z . The system can be 
characterised by the spin quantum numbers S and Ms- 

• For |<7i<7 2 ) we can use a product function. 
4.2.1 Electrons are indistinguishable 

In microscopic systems similar particles (e.g. electrons) can experimentally 
not be distinguished. The probability density ^/*^/ is in such many particle 
systems not allowed to be dependent on the arbitrary enumeration of the 
particles. A change of the numeration, exchange of particles, is not allowed 
to change 

Conclusion: With an exchange of two particles of the same type, the many- 
body wavefunction \l/ either stays the same or changes sign. Wave functions 
are either symmetric or anti-symmetric w.r.t. the exchange of two particles 
of the same kind. 

Spin eigenfunctions for two-electron systems: 

A naive product function leads to the following combination 

(a(l)a(2) 
a(l)/?(2) 
0(l)a(2) 
0(1)0(2) 

In the second and third product we can distinguish the two electrons, this is 
not possible for the products one and four. In the first case symmetric wave 
functions can be constructed through linear combinations: 

2- 1 / 2 [a(l)P(2) + P(l)a(2)] 
2- 1 /2[ a( l) /3( 2)- 0(1)^(2)] 
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An acceptable set of product functions is 



a(l)a(2); 

2-^[a(l)m - 0(l)a(2)]; 
/3(l)/?(2) 

Electron spin for two-electron systems: Singlets and triplets 



Spin eigenfunction 


Exchange 




5 


Notation 


a(l)a(2) 


symmetric 


1 


1 


triplet 


2- 1 /2 [a (i) /g( 2) + / g(l)«(2)] 


symmetric 





1 


triplet 


mm) 


symmetric 


-1 


1 


triplet 


2-^[a{l)m)- /?(l)a(2)] 


antisymmetric 








singlet 



Simple example: 

S z a(l)a(2) = (s zl + s z2 )a(l)a(2) 

= s zl a(l)a(2) + s z2 a(l)a(2) 

= -ha{l)a{2) + -ha{l)a{2) 
= ha{l)a{2) -> M 5 = 1 

4.3 Pauli principle 

All in nature realisable states of many-electron systems are described by 
space-spin-wave functions that are anti-symmetric w.r.t. electron exchange. 

Special case of the Pauli principle: Every spatial orbital can be occupied by 
at most two electrons with opposite spin. 

Generalisation to other particles: 

Fermions = Particles with half-integer spin 
Bosons = Particles with integer spin 



60 



All systems with several equal fermions (bosons) are described by wave func- 
tions, that are anti-symmetric (symmetric) w.r.t. exchange of two equal 
fermions (bosons). 

Example for fermions: electron, proton, neutron 
Example for bosons: H-Atom, He-nuclei 
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Lecture 5 

Hartree— Fock Approximation 



5.1 Slater-determinant 



We are looking for a n-electron wave function \1/ within the orbital model that 
satisfies the Pauli principle. This wave function is build from orthonormal 
spin orbitals 77. 

V(r,a) =tt(f)|<r) 



Short hand for electron 1: 

I7(l) = tt(l)|(7(l)> 

Ansatz 



tfa(l) 



*(1) 
*(1) 



1 



^l(l) 


^i(2) • 


• >?i(JV) 




%(2) • 


• V2 (N) 


^at(I) 


^(2) • 


■ Vn(N) 



Exchange of two electrons corresponds to an exchange of two columns. This 
leads to a change of sign for the determinant. The determinantal wave func- 
tions is therefore anti-symmetric w.r.t. an exchange of electrons (Pauli prin- 
ciple). 

A determinant vanishes if two rows are the same. In non-trivial wave func- 
tions \& 7^ all spin orbitals have to be different. A given spatial orbital can 
only appear twice in the determinant, combined with either a or /3. 
Pauli principle in the orbital model: A spatial orbital can at most be occupied 
by two electrons with opposite spin. 
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Writing out the determinant leads to N\ products of the form 771772 ■ ■ ■ i] N . All 
possible permutations of the N electrons in the N spin orbitals appear. In 
this approach the electrons are evenly spread out over all spin orbitals and 
therefore indistinguishable. 

Short notation: Often we omit the normalisation constant (AM) -1 / 2 and we 
only write the diagonal elements. 



*(l,2,...,N) = \r jl (l)r j2 (2)-.-r jN (N)\ 

Slater determinants are always eigenf unctions of S z , however, they are not 
necessarily eigenfunctions of S 2 . For the general case there are always linear 
combination of determinants that are eigenfunctions of S z and S 2 at the same 
time. Such spin-adapted linear combination of determinants (configurations) 
are needed to describe open-shell systems. 

Important special case: In closed-shell systems the spins of all electrons are 
paired. 



2,...,N) — 1^(1)^(2), tt 2 (3)* 2 (4), . . . * N/2 (N - 1)* N/2 (N) \ 



This Slater determinant is an eigenfunction of S z and S 2 with S = Ms = 
(Singlet). 



5.1.1 Energy expectation value for Slater determinants 

The electronic energy E e \ is the expectation value of the non-relativistic elec- 
tronic Hamilton operators H e \ w.r.t. the iV-electron wave function, written 
as a Slater determinant of orthonormal spin orbitals. 

E el = |# e l|*> 

N N 2 

i=l i<j l i 
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Using the Slater determinant function we get 



N N 



i=l i<j 



4° = / / ^(l)^.(l)£r ? ;(2)r ?l (2)rfr 1 rfr 2 



is the one-electron integral, J>™ Coulomb integral, and K^' an exchange 
integral. In these equations all integrals are over spin orbitals rj. By inte- 
gration over the spin coordinates we get to a description in terms of spatial 
orbitals. 

Example for the transition to spatial orbitals: 



A similar separation of spatial and spin parts occurs in all integrals. This is 
due to the fact that H c \ does not contain any spin coordinates. A systematic 
investigation of all integrals leads to: 




hu ] (a\a) = h\f ] 






r(*) 




equal 
opposite 



Interaction of two electrons in spatial orbitals ^j: 
Jij with opposite spin 

Jij — Kij with equal spin (smaller, because K i3 > 0) 



Important special case: closed shell systems 



N/2 



N/2 N/2 




i=l 



i=l j=l 
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5.1.2 Unitary transformation of spin orbitals 

In the definition of the Slater determinant we regard the spin orbitals rji(j) 
formally as elements Aji of a matrix A. 

= (N\y 1/2 det(A) 

A unitary transformation U of the spin orbitals leaves them orthonormal: 

i 

u- 1 = u\ (u- r ) lk = = {U*) ki 

In matrix form: 

A' = AU 

V = (Niy^detiA') 

= (N\)- 1/2 det(AU) 

= (N\)- 1/2 det(A) det{U) 
^' = ^ det(U) 

Calculation of det(U) for a unitary transformation: 

U ] U = 1 

det(U ] U) = det{U r ) det(fJ) = [det(U)\* det(U) = \det{U)\ 2 = 1 
det(U) = e l * 

^' differs from ^ at most by a phase factor e l * with ±1 for the real case. 
For all expectation values we get 

= (*'|6|*') 

In this sense determinantal wave functions are invariant w.r.t. unitary trans- 
formation of spin orbitals. 

5.2 Hartree— Fock method 

Goal: Calculation of optimal orbitals 

Approach: Variational method, minimisation of the energy expectation 
value for the Slater determinant under the condition that orbitals stay 
orthonormal. 
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Derivation for spin and spatial orbitals is similar. We study the more general 
case of spin orbitals. For orthonormal spin orbitals we have the energy in 
simplified notation: 

N N N 

E C 1 = ^2 h n + ^ E ~ K v) 

i=l i=l j=l 

Sij = (Vi\Vj) = &ij = $ij 



All integrals are for spin orbitals. Spin orbitals may be complex. Optimal 
spin orbitals minimise the energy expression while maintaining the constraint 
condition. This can be achieved by the method of Lagrange multipliers: 



N N 

C[{Vi}} = E el [{ Vi }} - E E e ^ - 6 «) 

i=l j=l 

5C = 



Lagrange multipliers have to be Hermitian in order for the Lagrange- 
Function C to be real. Variation of the spin orbitals must lead to a vanishing 
variation of C: 

N N 

5C = 5E cl € j^Sij = 

i=l j=l 

We investigate the variation of the integrals that is caused by the variation 
of the spin orbitals (c.c. = conjugate complex) 

SSa = j 5r 1 *{l)r lj {l)dT 1 + J i£(l)tfife(l)dn 

^ N N 

5E cl = E Shu + - E E^ _ 5K v) 

i i=l j=l 

6hu = J Srj;(l)h(l)rj j (l)dr 1 + J ^(1)^(1)^(1)^1 
8h u = J Sr]* (1)^1)^(^ + 0.0. 
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Because of Hermiticity of h, the second term in the equation for Shu is the 
complex conjugate of the first term. 



SJij = 



J J dr 1 dr 2 Sr ] *(l)r ]t (l)^-r ] *(2)r ]] (2) A 

J j ^^^(1)^(1)^(2)^(2) B 

J J dndr^* (1)^(1)^(2)^(2) C 

J J dr 1 dr 2 r ] *(l)r ]t (l)^-r ] *(2)Sr ]j (2) D 

. . A 1 A <f „ „ 2 

EE^ = EE 2 / /w^(i^(i)^(%(2)+c.c. 

1=1 7 = 1 1=1 7 = 1 J 12 



In the double sum in the last equation the terms A and C appear in equivalent 
pairs and can be combined (same is true for SK^): 



N N N N , , 2 

EE^' = EE 2 / / dr 1 dr 2 Sr ] ;(l)r ]3 (l)^r ] *(2)r ]l (2) + c.c. 

a — 1 a — 1 a—, J J 12 



i=l j = l i=l j = l 

For the remaining terms we get: 



N N N N 



E E e * 6S v = E E e ^ / h^hwt! + c.c. 



i=i j=i i=i j=i 



A further simplification can be achieved with the introduction of the Coulomb 
operator j, and exchange operator Kf 



3,(1)^(1) = j dr 2 r ] *(2)^r ]j (2)r ]i (l) 

J T\ 2 

K 3 (l) m (l) = [ dr 2V *(2)— 1^(2)^(1) 

J i3 = j ^1^(1)^(1)1^(1) 

Kij = J ^(1)^(1)^(1) 
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Taking all terms together we get: 



+ 



+ c.c. = 



3 3 

Variations Srj* are arbitrary and therefore every single integral to an index i 
has to disappear and the integrand [. . .] has to be zero. 



w + E^ 1 )-^ 1 )) 



V 



(1) = E^(1) fori = 1,..., N 



The same steps performed on the conjugate-complex part lead to the same 
result. 



Introduction of the Fock-operator: 

3=1 

F(1)^(1) = E^(1) 

3 

The current derivation did not lead to the desired eigenvalue equation for spin 
orbitals. We now investigate unitary transformations of the spin orbitals. 
We have already shown that determinantal wave functions and their expec- 
tation values are invariant w.r.t. such unitary transformations. The same is 
true for the Fock operator. For the Lagrange multipliers we get e^: 

(ifr(i)i£(i)iifc(i)> = E^<%(i)i^(i)> = c« 

3 

The multipliers are equivalent to the matrix elements of the Fock operator. 
A transformation of the spin orbitals leads to (in matrix form): 

e' = UhU 

e is an Hermitian matrix. Therefore there is a special unitary matrix U, 
which transforms e to diagonal form. For this special set of transformed spin 
orbitals we get 

*"(l)ifc(l)=e;ifc(l) 

Normally one omits the special mark for these canonic spin orbitals and one 
writes: 

Fr]i = eir]i for i = 1, . . . , TV 
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5.2.1 Solution of the Hartree-Fock equations 

Problem: The Fock operator consists of Coulomb Jj and exchange operators 
Kj. Therefore, it depends on the spin orbitals rji, we want to calculate 
with the Hartree-Fock equations. 

Solution: Iterative method starting from an arbitrary but reasonable initial 
set of orbitals r]f\ 

Hartree-Fock equations have to be solved iteratively. 
Scheme 

R (0) } — FW — — _> _> etc. 

This sequence is stopped, if the solutions in two consecutive iterations don't 
change (within a given tolerance). This method is called SCF (Self-Consistent- 
Field). 

For a known (self-consistent) set of occupied spin orbitals the Fock operator F 
is a well-defined Hermitian operator with an infinite number of eigenfunctions 
rji (N occupied spin orbitals with lowest orbital energies and un-occupied 
(virtual) orbitals with higher orbital energies). 

5.2.2 Orbital energies and Koopmans' theorem 

For an iV-electron system, minimisation of the energy of the determinant 
= 1 771772 • • 'Vn) leads to an eigenvalue equation J-\r] a ) = e a \i] a ) for the n 
occupied spin orbitals {tja}- The Fock operator has a functional dependence 
on these occupied spin orbitals, but once the occupied spin orbitals are known 
the Fock operator becomes a well-defined Hermitian operator, which will have 
an infinite number of eigenfunctions 

= £j\Vj) J = 1,2, ...,00 

Each of the solutions has a spin orbital energy e r The N spin orbitals with 
the lowest orbital energies are just the spin orbitals occupied in |\&) for which 
we use the indices a,b, . . . . the remaining infinite number of spin orbitals 
with higher energies are the virtual or un-occupied spin orbitals, which we 
label with the indices r, s, . . . . 

The matrix elements of the Fock operator in the basis of the spin orbital 
eigenfunctions is diagonal with diagonal elements equal to the orbital energies 

(Vi\F\ Vj) = ej{Vi \Vj) = tjtij 
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The orbital energies can be expressed 



= (i\h\i) + 66) - (i6, , 

6 

where we have introduced the short form (ij, fcZ) for the two-electron integrals 

(ij,kl) = J ^^^(1)^(1)^(2)17,(2) . 

In particular we get for occupied orbitals 

e ffi = (a | h | a) + y^((aa, 66) — (a6, 6a)) 

= (a | a | a) + y^((aa, 66) — (a6, 6a)) 

and for virtual orbitals 

e r = (r \ h \ r) + ^~]((rr, 66) — (r6, 6r)) 

The orbital energy e a represents the energy of an electron in the spin orbital 
1 77 a ). This energy is the kinetic energy and the attraction to the nuclei ((a | 
h | a)) plus a Coulomb and exchange interaction with each of the remaining 
N — 1 electrons. The result for e r has a different character. It includes the 
kinetic energy and the nuclear attraction of an electron in \r] r ), but includes 
Coulomb and exchange interactions with n electrons of the Hartree-Fock 
ground state. It is as if an electron has been added to |\& ) to produce a 
(N + l)-electron state. 

If we simply add up the orbital energies of the occupied states, we get 

N N N N 

a a a b 

If we compare this with the total Hartree-Fock energy 

N ^ N N 

E HF = ^2 (a | h | a) + - ^ ^ ((aa, 66) - (ab, ba)) 

a a b 
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we see that 



JV 



EKF 7^J2 € ' 



a 



a 



and the total energy of the state |^o) is not just the sum of the orbital 
energies. The reason is as follows. The energy e a includes Coulomb and 
exchange interactions between an electron in orbital a and electrons in all 
other occupied spin orbitals, in particular orbital b. But includes Coulomb 
and exchange interactions between an electron in orbital b and electrons in 
all other occupied spin orbitals, in particular orbital a. Thus if we add e a 
and €b we include the electron-electron interactions between an electron in a 
and one in b twice. The sum of orbital energies counts the electron-electron 
interactions twice. This is the reason for the factor | in the correct energy 
expression for the total energy relative to the sum of orbital energies. 

If the total energy is not the sum of the orbital energies what physical sig- 
nificance can we attach to orbital energies? 

Koopmans' theorem 

Given an iV-electron Hartree-Fock single determinant l^^o) with occupied 
and virtual spin orbital energies e a and e r , then the ionisation potential to 
produce an (N — l)-electron single determinant | Ar ~ 1, I' a ) with identical spin 
orbitals, obtained by removing an electron from spin orbital f] a , and the 
electron affinity to produce an (N + l)-electron single determinant \ N+l ty r ) 
with identical spin orbitals, obtained by adding an electron to spin orbital 
r] r , are just — e a and — e r , respectively. 

Proof: The ionisation potential is 



where N 1 E a and N E Q are the expectation values of the energy of the two 
relevant single determinants 



IP = 



JV-l 








N 



E a = ( N ~^ a i n r- 1 



JV-l 
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The ionisation potential for orbital c is therefore 
IP = N ^ E c - N E 

1 X 

= ^(a | h | a) + - ^ XI ^ aa ' && ) ~~ ( a& ' &a )) 



2 

- J^(a | /i | a) + ^ X X ((aa, 66) - (ab, ba)) 

a a b 

= -(c | h | c) - - X ((aa, cc) - (ac, ca)) — - ((cc, 66) - (6c, cb)) 

a b 

= -(c\h\c}-J2 (( cc > bb ) ~ ( fcc - cb )) 

b 

= -€ c 

Occupied orbital energies are generally negative and ionisation potentials are 
positive. 

With a very similar calculation we can also calculate 

EA = N E - N+1 E r = -e r 

This result is consistent with the previous observation that e r included inter- 
actions with all N other electrons of the ground state and thus describes an 
N + 1th electron. 

Koopmans' ionisation potentials are reasonable first approximations to ex- 
perimental ionisation potentials (missing are orbital relaxation and corre- 
lation effects). Koopmans' electron affinities are unfortunately often bad. 
Many neutral molecules will add an electron to form a stable negative ion. 
Hartree-Fock calculations on neutral molecules, however, almost always give 
positive orbital energies for all the virtual orbitals. 

5.2.3 Brillouin's theorem 

The Hartree-Fock equation produces a set of spin orbitals. The single de- 
terminant |^o)) formed from the N lowest spin orbitals is the Hartree-Fock 
approximation to the ground state. There are many other determinants with 
N functions that can be formed from the full set of spin orbitals. Brillouin's 
theorem makes a statement of a subset of these determinants. This subset 
is the set of singly excited determinants \^/ r a ) obtained from \^/ ) by a single 
replacement of r\ a with r\ r . 

Brillouin's theorem 

Singly excited determinants will not interact directly with a reference 
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Hartree-Fock determinant |^ ), i-e., (^o I H I = . 
Calculating the matrix element 

(tf o I U | %) = (a | /i | r) + ^ ((or, 66) - (aft, 6r)) 

The right hand side of this equation is just the matrix element of the Fock 
operator in the spin orbital basis. Therefore, 

(tf I H | K) = (Va I T | Vr) 

The matrix element that mixes singly excited determinants with \ty ) is thus 
equal to an off-diagonal element of the Fock matrix. Now, by definition, solv- 
ing the Hartree-Fock eigenvalue problem requires the off-diagonal elements 
to satisfy (rji \ T \ rjj) — (i ^ j). One can then say that solving the 
Hartree-Fock eigenvalue equation is equivalent to ensuring that |^o) wm n °t 
mix with any singly excited determinants. 
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Lecture 6 

Molecular Orbital Theory 



6.1 LCAO-Ansatz 

LCAO: Linear Combination of Atomic Orbitals 

In the LCAO approach the molecular orbitals (MOs) ^ j are written as linear 
combinations of atomic orbitals (AOs) 

In the following the equations can also be applied to other types of basis sets 

that are not necessarily of atomic type. 

Conventions: 

k,l, . . . for MOs 
fi, v, X, a, . . . for AOs 

We restrict the following derivations to closed-shell systems, singlet states 
with S = M s = 0, with all orbitals doubly occupied. 
Energy formula: 

N/2 N/2 N/2 

i=i i=i j=i 

K = j ^(1)^(1)^(1)^ = (VilhVi) 

Jij — J y'*r(l)* i (l)^(2)* i (2)dr 1 dr 2 = (w,^) 
K%j = (U>U) 

We insert the LCAO equation and express all MO integrals with AO integrals. 
We assume that MOs and AOs are real functions. 
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Connection of MO and AO integrals: 

h ii = (y~] c ^^v i h E Cyi ^ v ) = E E c ^°^ h 



fl V I /J, V 



Ji 'j = E E E C ^ C ^ C ^ C Tj(^ Af7 ) 

|U f A CT 

^ = E E E E fyWKC'jiM X(J ) 

jl V \ o 

Xa) = J J $ M (l)$ i ,(l)^$ A (2)$ CT (2)rfr 1 rfr 2 
Energy expression in AO basis: 

-Eel = 2 ^ ^ ^ ^ ^ ^ c pii c uih^u 
i fi v 

+ E E E E E E c ^ c ^CAjC (Ti [2(^Z/, Act) - (//A, z/<r)] 

i j H v \ <t 



Introduction of the density matrix P^ u : 



N/2 



Simplification of the energy expression in the AO basis: 

E * = E E p ^ v + ^ E E E E ^^[2(^, A(t) - (//A, zat)] 

For a given basis all the AO integrals and (fiu, Act) are (in principle) 
known. 

The LCAO coefficients c M « and the density matrix elements P^ are not 
known. They can be calculated from a linear variational calculation (Root- 
haan-Hall method). 



6.2 Roothaan-Hall equations 

Instead of a complete derivation it is easier to insert the LCAO approach 
in the general Hartree-Fock equations and multiply from left with $* and 
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integrate: 



Introduction of the Fock matrix F and the overlap matrix S in AO basis: 

F, u = (%\F$ U ) 
S, u = 

Eigenvalue equation: 

^ ^ FfjiuCut £j ^ ^ S^ u c u i i 1,2,..., 

FC = SCe Matrix-Notation 

F = {F^}, C = {c vi }, S = {S^}, e = {eAj} 

For K basis functions $ M there are K equations. The matrices have the 
dimension K x K; F and S are symmetric, e is diagonal. 
Fock matrix elements in AO basis: 

N/2 

F, v = ^,\h + Y J ^J ] -K j )\^> u ) 

i=i 

N/2 

i=i 

JV/2 

j=l A a 

F^ = + X X p A CT [(/^, A(J ) ~ ^(/ iA > az/ )] 

A cr 

Solution of the Roothaan equations 

As for the general Hartree-Fock equations, the Roothaan equations can only 
be solved iteratively: The Fock matrix elements F M „ depend on the density 
matrix elements P\ a . The density matrix has to be calculated from the 
unknown LCAO coefficients c„j. 
A LCAO-MO-SCF calculation: 

1. Choose a basis of atomic orbitals {$ M } 
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2. Calculate the AO integrals h^ u , (/lis, Act) and S^ v . 

3. Choose a reasonable initial density matrix {Pfk)}. 

4. Calculate the corresponding Fock matrix. 

5. Solve the eigenvalue equation (Roothaan-Hall). 

6. Calculate a new density matrix. 

7. Iterate (4)- (5)- (6) until self consistency is achieved. 
Solving the eigenvalue equations: 

We examine a linear transformation of the non-orthogonal basis {^j into 
an orthogonal basis 

V 

A a 

= E E X X^XaX av 
A a 
I = X ] SX 

The overlap matrix S in the basis is Hermitian and can be diagonalised 
using a unitary transformation U . 

U^SU = s 

The diagonal matrix s contains the eigenvalues of S (sa > 0). For the 
transformation matrix X the following approaches are used: 

X = S^ 1 ^ 2 = Us~ l l 2 U^ symmetric orthogonalisation 

X = Us^ 1 / 2 canonical orthogonalisation 

The symmetric orthogonalisation is used most often. The canonical orthog- 
onalisation has the advantage that linear dependencies in the basis (sa ~ 0) 
can easily be eliminated. 

Connection between the different sets of LCAO coefficients: 



= E = E ii E x ^ = E 

c = xc', d = x~ x c 



Xu ^°^i 
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With these definitions we can rewrite the eigenvalue equations: 

FXC' = SXC'e 

(X^FX)C' = (X^SX)C'e = C'e 

Using the Fock matrix F' in the basis {$^} we get a standard eigenvalue 
equation: 

F' = X^FX 
F'C' = C'e 

Modified LCAO-MO-SCF procedure: 

1. Choose a basis of atomic orbitals {$^} 

2. Calculate AO integrals h^, {jiu, Aer) and S^. In addition calculate the 
matrix X. 

3. Choose a reasonable initial density matrix {P$}. 

4. Calculate the Fock matrix. 

5. Solve the Roothaan-Hall equations. 

(a) Transform the Fock matrix 

(b) Solve eigenvalue equation 

(c) Back transform coefficients 

6. Calculate a new density matrix 

7. Iterate (4)-(5)-(6) until self-consistency is reached 

6.3 Unrestricted open-shell Hartree— Fock 

Up to now we have first derived the Hartree-Fock equation for general spin 
orbitals rj 

and then have had a look at the special case of closed shell orbitals. The 
special case lead to the restricted Hartree-Fock (RHF) equations. This spe- 
cial case has the advantage that the wavefunction is an eigenfunction not 
only of the Hartree-Fock Hamiltonian, but also of the spin operators S z and 
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S 2 . For the case of unpaired spins there are two possible methods avail- 
able. Either the electrons are partitioned in paired and non-paired sets and 
the paired electrons treated with equal spatial orbitals (the restricted open- 
shell Hartree-Fock approach, ROHF) or all electrons are treated on the same 
footing using general spin orbitals in an unrestricted Hartree-Fock approach 
(UHF). The first method has the advantage that the wavefunction is still 
an eigenfunction of the spin operators but the disadvantage that it needs a 
rather involved theory and that the solution is not variationally stable w.r.t. 
an unrestricted wavefunction. The UHF method is conceptually simpler and 
leads to the variational ground state, but might lead to wavefunctions that 
are no longer eigenfunctions of S 2 . 



6.3.1 Spin contamination 

We start from a Slater determinant 

I^uhf) = \V1V2 ■•■Vn) 

where the spin orbitals are 



??(x) 



either occupied with an a or j3 electron. There are N a a and Np f3 electrons 
in the system. 

N = N a + N p 

If we assume that N a > Np then we can easily find that the expectation 
value of S z is 

= \ (N a - N P ) 

The exact expectation value of a corresponding wavefunction for S 2 would 
be 

(SLct) = \ (N a - N p ) (N a -Np + 2) 

However, if we calculate the expectation value for an unrestricted Slater 
determinant, we get 

(^H F ) = (^xact)+^-EEl^ /3 | 2 

i j 

where 

sf = <C I 

The value for (<Suhf) should always be monitored when working with unre- 
stricted wavefunctions. 
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6.3.2 Pople— Nesbet equations 

The Hartree-Fock equations for a spin orbital is 

where a stands for a general spin function. If we multiply from the left with 
a(cui) and integrate over the spin variable uj 1 we get 



^ = M 
where the spatial Fock operators are defined by 

^"(rx) = J du 1 a*(uj 1 )J : '(r 1 ,uj 1 )a(uj 1 ) 

The operator jF a (l) is the kinetic energy, nuclear attraction, and effective 
potential of an electron of a spin. The effective interactions of an electron of 
a spin include a Coulomb and exchange interaction with all other electrons 
of a spin plus only a Coulomb interaction with electrons of j3 spin. Thus 

= h(i) + e [^(i) - ^(i)] + e Jsf(i) 

a a 

The sum over the N a orbitals formally includes the interaction of an a elec- 
tron with itself. However, since 

[J?(l) - J^(l)] #(1) = 

this self-interaction is eliminated. The corresponding Fock operator for elec- 
trons of (5 spin is 

= mi) + E W 1 ) - kP M + E J a (!) 

a a 

The kinetic energy and nuclear attraction of an electron in one of the unre- 
stricted orbitals <pf or 0f is the expectation value 

h« = (€\h\<P?) or hi = {tf\h\tf) 
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The Coulomb interaction of an electron in (j)f with one of opposite spin in 
0f is 

Jtf = jfr = m i j? i #> = i j: i ^) = («, 

The corresponding Coulomb interactions between electrons of the same spin 
are 





c i ^ 


C> = < 








1 4) = < 


# 1 4 1 



and 

The exchange interactions between electrons of parallel spin are 
and 

K?f = (0f | < | 4) = I K* | 0f ) = (0f , 0f0f ) 

There is no exchange interaction between electrons of opposite spin. 
The total energy for a unrestricted determinantal wavefunction is 



a a 



a 6 



^,9 W/3 V a ^ 

a/3 
ab 

a b 



The factor of 1/2 in the third and fourth terms eliminates the double counting 
in the free summation. The self-interaction disappears since J" a Q — = 

Basis set expansion of the orbitals (LCAO) 

K 

€ = Y, C ^ i = l,2,...,K 

K 

^=eo. <=i,2,...,^ 

Notice that both sets of orbitals are expanded in only one set of basis function. 
The two eigenvalue equations for the orbital sets guarantee that the eigen- 
functions for a and (3 sets individually form orthogonal sets. There is no 
reason, however, that a member of one set need be orthogonal to a member 
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of the other set. Even though the two sets of spatial orbitals overlap with 
each other, the set of 2K spin orbitals will form an orthonormal set, either 
from spatial orthogonality (act and (3(3 case) or spin orthogonality (a(3 case). 
Substituting the expansions of the orbitals into the Hartree-Fock equations 
gives 

V V 

If we multiply this equation by Xu(l) an d integrate we get 

E F »v C »i = e °i E S i^ C lj 

V V 

where S is the overlap matrix of the basis functions and F CT is the matrix 
representation of T° 

f; v = J dv lX i(i)^{i) Xu {\) 

Therefore we end up with two coupled matrix equations, the Pople-Nesbet 
equations 

F a c a = SC « e a 



Unrestricted density matrices 

The charge density of the system is 

p T « = E i^wi 2 + E i^«i 2 = p a ^ + ^« 

a=l a=l 

The corresponding spin density is defined by 

/(r) = p«(r)-/(r) 
By substituting the basis set expansion for the orbitals we get 

P a (r) = ElCWI 2 = EE P ^Wx,(r) 

<2=1 /I V 

and 

Ng 

A*) = Ei^wi 2 = EE p ^«^ r ) 

a=l /j, v 
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where the a and (3 electron density matrices are used 

pa \ ^ /~ict I /~ia \* 

a 

N/3 

p/3 _ \^ f(~if3 \* 
fiu / j [ia\^ ua) 

a 

In addition to these two density matrices we can also define a total and 
density matrix 

pT _ pa _|_ p/3 
-pS _ pa p/3 

With these definitions we can write the Fock matrices in compact form 

= k »T + E E ffA ) " P A.(M, OV) 

A cr 

F i> = KZ* + E E P ^ (7A ) " P A.(M, (71/) 
A cr 

and the total energy is 
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Lecture 7 



Correlation Energy 



7.1 Basics 



For many systems the Hartree-Fock method accounts for about 99 % of the 
total energy. However, often the remaining 1 % are important for physical 
and chemical properies of the system. 
Correlation: 

Reduced probability to find an electron with different spin close to another 
electron: Coulomb hole. 

Reduced probability to find an electron with the same spin close to another 
electron: Fermi hole. 

Hartree-Fock: Variational optimal one determinant wave function. 
Improvements are only possible with multi determinant wave functions. 



Produce new determinants by exchange of MOs in $hf with one, two, three, 
. . . virtual MOs. Singly, doubly, triply, . . . excited determinants. 




Excited Slater determinants: 
Hartree-Fock determinant 



• • • $JV/2$iV/2 



Not used MOs (virtual, unoccupied MOs) 



^ N/2+l, • • • , §K 



85 



This procedure generates determinants with all possible eigenvalues of S z , 
linear combination of these determinants that are also eigenfunctions of S 2 
are called CSF (Configuration State Functions). 
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The complete optimization problem to find an exact wave function with this 
method is two dimensional. Enlargement of the one-electron basis to get 
a converged Hartree-Fock wave function also increases the number of vir- 
tual MOs. This means an increased number of excited determinants. The 
complete solution for a given basis set (with K functions) is called Full Con- 
figuration Interaction (FCI). 



Basis 
complete 



-+ 



minimal 



exact HF 



exact 
— |— solution 



HF 



FCI 
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7.2 Matrixelements 

Examples of determinants with two orbitals 

% : -L [a(l)6(2) - 6(l)a(2)] 

^: ^ [r(l)6(2) - 6(l)r(2)] 

Kt ■■ ^ [r(l)a(2) - s(l)r(2)] 



($o | O | $ ) = ^<a(l)6(2) - 6(l)a(2) I O | a(l)6(2) - 6(l)a(2)> 

= \{a{l)b{2) | 6 | a(l)6(2)> - \(b{l)a{2) \ 6 | 6(l)a(2)> 
- J(a(l)6(2) | 6 | &(l)a(2)) + J<6(l)a(2) | 6 | a(l)6(2)) 



Operators: 



Case 1: 



Case 2: 



O = 1 

AT 



i=i 

N N n 

e z 



($1$) = i 

(^|O2|$> = ^E(i<Uj)-(v'U0 



($|$;> = o 

($|Oi|$;> = (a\h\r) 
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Case 3: 

(*|di|*2> = o 

($|6 2 |« = (ar\bs) - (as\br) 

Case 4: 

(*|0 2 Kc) = o 

Energies of determinants: 

i 

/iji : contribution from each electron 

Ej : contribution from each unique electron pair 

Ek '■ contribution from each unique electron pair with same spin 

Examples 

Orbital 

2 +- f- - 44 

■ f 4 # f 

abed 

a) /in + h 2 2 + J12 - Kn 

b) h u + h 2 2 + J 12 

c) 2h u + Jn 

d) h u + 2h 22 + 2J 12 + J 22 - K 12 

7.3 Configuration interaction method (CI) 

Ansatz 

*ci = a $HF + a s® s + Yl aD ^ D + " ' 

S D 

Variational calculation with normalisation constraint 

(^cil^ci) = 1 
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Lagrange function 

C = (tt C i|£|*ci) - A [(*ci|*ci) - 1] 
(*ci|ff|*ci> = X;E a < a ^ $ <l^l^> 

j=0 j=0 

= + y ^y£ d a i a j ^ i \H\^ j ) 

i=0 i=0 jVj 

(*ci|*ci) = X)I] a<a ^ $i l^) 

i=0 j=0 
i=0 

=5>. 2 

i=0 



Variation: 



<9£ 



— = 2^a J ($ l |i7|$ J )-2Aa J = 

a t |^> - A) + ^a^if |$,-> = 

a, (Ei- A) + J2 = 

Multiplication of the last equation with a, and summation over all equations 
leads to 

i 11 jy^l 

and therefore 

A = (*ci|#|*ci) • 

The Lagrange multiplier is the CI energy (A = E C i). We have a CI eigenvalue 
equation in matrix form 

Ha = £a 

Size of the CI matrix: 

N electorns (indistinguishable) in K MOs. 

, K \ K\ 

Number of Determinants = 
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N J N\(K-N)\ 



Many determinants have a wrong (not desired) M s and can be neglected. In 
addition the CI matrix is reduced in the CSF basis 

K\(K + 1)! 

Number of CSF(Singlet) - 



(N/2)\(N/2 + 1)\(K - N/2)\(K - N/2 + 1)! 



Examplel: 

H 2 10 electrons in 38 AOs (= 38 MOs) 



Excitation CSF total CSF 






1 


1 


1 


71 


72 


2 


2485 


2557 


3 


40040 


42597 


4 


348530 


391127 


5 


1723540 


2114667 


6 


5033210 


7147877 


7 


8688080 


15836557 


8 


8653645 


24490202 


9 


4554550 


29044752 


10 


1002001 


30046752 



Sizes: 

Cl-vector w 30 • 10 6 elements (« 240 Mbytes) 

Cl-matrix « (30 • 10 6 ) 2 = 900 • 10 12 elements (« 7200 Tbytes) 

but most Hij = 0, number of non-zero matrix elements is proportional to the 

matrix size. 

Solution of the CI eigenvalue equation: 
Only lowest eigenvalues are of interest. 
Iterative diagonalization for some of the lowest eigenvalues 

• Lanczos method 

• Davidson method 

• Residuum methode 

The CI matrix is sparely occupied. The number of matrix elements ^ 
is proportional to the dimension of the matrix. H is not explicitely calculated 
but only Ha. 
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Reduced CI methods: 

Idea: Not all determinants are equally important. 

Ansatz: Only allow excitations from a subset of MOs into a subset of virtual 
MOs (active space). 

Allow only a maximal number of excitations. 



Ne-Atom: Weight of CSF 



Excitation 


Weight 





0.9644945573 


1 


0.0009864929 


2 


0.0336865893 


3 


0.0003662339 


4 


0.0004517826 


5 


0.0000185090 


6 


0.0000017447 


7 


0.0000001393 


8 


0.0000000011 


Active Space: restriction of orbitals for excitations 



active 
space 
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7.4 Dynamical and non- dynamical electron 
correlation 



For the majority of molecules, for example all moleculaes that can be assigned 
a single Lewis structure, the main error from in the Hartree-Fock approxi- 
mation comes from ignoring the correlated motion of each electron with all 
the other electrons. This is called dynamical correlation because it refers to 
the dynamical character of the electron-electron interactions. Even though 
the correlation energy may be large, the Hartree-Fock determinant domi- 
nates the configuration interaction expansion in these cases. The correction 
is made up from many small contributions from excited state determinants. 
In certain cases, however, it happens that a set of nearly degenerate orbitals 
occure in a molecule that can only be partally occupied. In a restricted 
Hartree-Fock calculation one or several of these orbitals cannot be included 
into the wavefunction. An equally good wavefunction could be generated by 
including these orbitals and exculding the others. In a subsequent config- 
uration interaction calculation we expect that these determinants will have 
similar weights. The correlation energy gained for these cases is called non- 
dynamical correlation. 
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Lecture 8 



Coupled Cluster Approaches 

8.1 Size Consistency 




Method X is size consistent if 

E x (R = oo) = 2E x (F) 

Example: Hartree-Fock 

^ r= oo = V Fl ■ ^F 2 

therefore E(R = oo) = E{F 1 ) + E(F 2 ) = 2E(F) 

Example: SDCI 
Minimal basis set (#2)2 
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R 



a 



Determinants for the full calculation: 



HF 


aabb 


singly 


arbb 


singly 


aabf 


CT 


asbb 


CT 


aabf 


doubly 


rfbb 


doubly 


aass 


CT 


asbf 



For R = oo the CT determinants are zero and we are left 



with 5 determinants. 

For the calculation of the separated systems we have three determinants per 
system. 

HF aa 

singly af 

doubly tt 

TT _. , 7 The total wavefunction will be a linear combination of all the 
tit bb 

singly bs 

doubly ss 

products. We get 9 determinants, including triply and quadruply excited 
ones. The energy of these two calculations will certainly differ. Truncated 
CI methods are not size consistent. 

2E(F) < E(R = oo) 

The bigger the system the smaller the percentage of correlation energy is 
included in a truncated CI calculation. 

E CoiI (SDCI)/E Cor r(FCI)^0 for N el -> oo 

A related concept is size extensivity. A method is size extensive if 

Ecoii(Fn) OC N ■ Ecorr(F) 
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8.2 Coupled Cluster Methods 

Perturbation theory based methods add all corrections (singly, doubly, triply 
. . . excitations) up to a given order to the reference function. In Coupled 
Cluster (CC) methods for a given type of correction all order will be included. 
Ansatz for the CC wave function 

^cc = e T $ 

oo 1 

e f = V if k 
^ k\ 

k=0 

The operator T is defined as 

f = f 1 + f 2 + f 3 + --- + f N , 

where the operator Tj applied to the Hartree-Fock wave function generates 
all i-th excited Slater determinants. 

occ vir 

a r 
occ vir 

a<b r<s 

The expansion coefficients t are called amplitudes. They are our new varia- 
tional parameters. 

e f = i + f 1 + (f 2 + \fD + (f 3 + f 2 f ! + If*) 

2 o 

+ (f 4 + f 3 fi + \fl + if 2 f f + if 4 4 ) + • • • 

Insertion into the Schrodinger equation: 

He f % = Ee f % 

multiplication from left with $g and integration 

(<& \He T %) = E cc (®o\e T <Z>o) 

= E cc {Q \(l + f 1 + f 2 + ---)Q ) 

= -^cc 
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We make use of the fact that the Hamilton operator only contains one- and 
two-electron terms: 

E cc = ($o|# (1 + Ti + f 2 + ^f?)® ) 
= ($o\H$ ) + ($ |#Ti$ > 

+ ($o|^T 2 $ ) + i($ |^f 1 2 $ ) 

occ vir 



a r 
occ vir 

For $ = $hf we have ($o|-^^a) = (Brillouins theorem) and the second 
integral can easily be written with MOs: 

occ vir 

E C c = E + J2J2 (*S + #2 - OS) [KM - (as\br)} 

a<b r<s 

The CC energy is completely determined by the Hartree-Fock energy and 
the one and two electron amplitudes. 

Coupled (nonlinear) equations for the amplitude can be generated from the 
Schrodinger equation by multiplication from left with $^*, . . . and sub- 
sequent integration. 

($ r a \ — > N occ ■ N virt equations 

Ml — ► N 2 occ ■ N 2 irt equations 

8.3 Approximate Coupled Cluster methods 

For the exact CC method all operators Ti{i < N) have to be used. 

f = Ti + f 2 + f 3 + • • • f N 

The method scales like M 2N+2 for M basis functions and iV electrons. The 
only practical method is to use approximations of the form T ps T\ + T 2 
(CCSD). 

f = f 1 + f 2 CCSD M 6 

f = f 1 + f 2 + f 3 CCSDT M 8 

f = fi + f 2 + f 3 + f 4 CCSDTQ M 10 
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The case T — Ti does not improve on the Hartree-Fock energy. Multiplica- 
tion of the Schrodinger equation from left with $q, and integration 
leads for CCSD to a system of coupled, non-linear equations for the ampli- 
tudes t r a und t^l- These equations have to be solved iteratively. The most 
common improvement over the CCSD method is achieved by a perturba- 
tive correction for T 3 . Inclusion of T 3 with perturbation theory leads to the 
CCSD(T) method with M 7 scaling. 



+ [fA + hrf^j + Qf| + l -f 2 fl + 1^ + • • • 



CCSD equations 



($ r a I if (l + fi + (V 2 + + ^f 2 fi + ^f^j $ > = #CCSD<^ I Ti$ > 

($3 I ff (•••) $ o) = e C csd(KI I (r 2 + \fA $ ) 



CCSD is non- variational. CCSD is size-consistent; higher excitations are 
included through products (t^f^). 
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Table 8.1: H 2 at the basis set limit 

Exp. CCSD(T) HF 



r OH [pm] 
a(HOH) 

fie [D] 

Atomization energy [kJ/mol] 



95.72 


95.64 


94.0 


104.52 


104.2 


106.2 


1.8473 


1.853 


1.936 


975.3 


975.5 


652.3 



Table 8.2: Basis set error: Statistics over many small molecules 



Basis 


HF 


CCSD - HF 


CCSD(T) - 


cc-pVDZ 


5.8 


23.8 


1.3 


cc-pVTZ 


3.9 


6.3 


0.5 


cc-pVQZ 


1.5 


1.5 


0.2 


cc-pV5Z 


0.2 


0.5 


0.0 


cc-pV6Z 


0.0 


0.4 


0.0 



.4 Accuracy in Coupled Cluster calculations 

property 



A 
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Lecture 9 



M0ller— Plesset Perturbation 
Theory 

9.1 Electron correlation with perturbation 
theory 

Problem: What is and what is H {1) ? 
Perturbation theory 




#(0)^(0) = £(0)^(0) 
n n n 



First order energy correction to ground state ty g 



Second order energy correction 



and first order correction to wavefunction 



*?' = E 




and we have 



4 2> = E 
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with 



Definition of Hamiltonian for the unperturbed system 



ff (0) = E^) 



The unperturbed Hamiltonian is chosen to be the sum of all Fock-operators 
of all electrons of the system. This is a sum of one-electron operators with- 
out coupling. Therefore the solution will be the sum of all energies of the 
independent systems and the wavefunction will be the product of the one- 
electron wavefunctions (orbitals). Because of the Pauli-principle we have to 
anti-symmetrize this product and get a determinantal wavefunction. 



The complete basis of all zeroth-order solution are all possible determinants 
generated from the orbitals. 



The energy of the other determinants (basis functions) is also just the sum 
of the orbital energies. 



¥^ = |$ 1 $ 2 ---$;v| = $HF 



^ = |$ 1 $ 2 ---$ r ---$ JV 



and all other excited determinants ® r a f c , ■ ■ •)• 
Energy 



4 0) = E e ^^ 




= ei + e 2 + e 3 H 



Perturbation Hamiltonian 




h - #(o) = j2 h(i) + v cc -J2 P ^ 



i i 
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M0ller-Plesset perturbation theory 

MPO = £(MP0) = ($ H f|# (0) |$hf> = X) e< 

i 

MP1 = E(MPO) + E(MP1) = ^ ei + ($ hf |^ (1) |$hf) 
1 e 2 



«7 



i ij ij 



n 



= E(HF) 
MP2 = E(BF) + £(MP2) 

From general perturbation theory: 

E( MP2) =£ (^yw 

With 

($o|i/ (1) |$^) = (Brillouin Theorem) 
(%\H^\^ r a s b ) = (ab\rs) - (ar\bs) 
($ |^ 1 )|$^ c ) = 
and the energy denominator 

E - E($ r a l) = €a + e b -e r -e s 

we get 

E(MP2) = ££ [Hrs) - (ar|6s)]2 



a<b r<s 



MP perturbation theory is size-consistent and size-extensive. 

Order Scaling % correlation 



1 


M 4 





2 


M 5 


« 90 


3 


M 6 


« 95 


4 


M 7 


« 98 


5 


M 8 




6 


M 9 
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Table 9.1: Water, DZ basis set, RHF calculation 



Method 


r e 


1.5r e 


1.5r e 


MP2 


92.9 


86.9 


80.6 


MP3 


95.5 


88.0 


78.5 


MP4 


99.2 


96.8 


94.0 


MP5 


99.6 


97.5 


95.4 


MP6 


99.9 


99.3 


99.6 


MP7 


99.95 


99.5 


99.9 


MP8 


99.98 


99.9 


101.3 


FCI 


-0.130085 


-0.191996 


-0.299864 



MP2 efficient (~ HF) implementation possible 
recovers xs 90 % of correlation energy 

MP3 often not more accurate than MP2 
— > not a standard method 

MP4 improves on MP2 

similar cost, but inferior to CCSD(T) 

MP5- • • only for benchmarking 
Only MP2 is a commonly used method in quantum chemistry. 



9.2 Convergence of MP series 

No smooth convergence, or no convergence at all is possible. 

Detection of difficult cases: HF determinant is not a good starting point. 
Weight of $ in perturbation wavefunction is less than 90 %. Example: H 2 
at a 2.0r e geometry. The HF determinant has only a 58% weight in MP2. 
The wavefunction has multi-configuration character. 

MP based on unrestricted Hartree-Fock 

• Slow convergence due to spin contamination problems 

• To be used with care! Better use spin-projection methods. 
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MP based on multi-configurational wavefunctions (e.g. CAS-PT2) 



• expensive 

• often very accurate 

• needs expert knowledge 

9.3 The dissociation problem 

Restricted vs. unrestricted solutions. 

Example: H 2 — > H + H, homolytic bond breaking 
Minimal basis $1, $ 2 with overlap S 12 = ($1 | $ 2 ) 



Symmetry therefore fully defines the restricted solution \I>o = l^i^i 
Unrestricted case 



= [2(1 + s 12 )r i/2 ($! + $ 2 ) 

* 2 = [2(1 - ,S 12 )]- 1/2 ($! - $ 2 ) 



cos 6^1 + sin 6 * 2 
cos 0*i- sin 9 * 2 



if 6 = then RHF=UHF. 
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Bond breaking (homolytic) 

• restricted wavef unctions fail 

• unrestricted wavefunctions may have spin-contamination 

• MP2 is highly affected 

• CCSD(T) might solve problems partially 
Ultimate solution: Multi-configurational SCF (MC-SCF) 

|^o) = a $o + + • • • 

The wavefunction is approximated by a small CI expansion. Variationally 
optimize orbitals and Cl-coefficients at the same time. 
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Lecture 10 



Density Functional Theory: 
Part I 

10.1 Basics 

In the following we investigate iV-electron systems within the Born-Oppenheimer 
approximation. The system is described by a Hamilton operator H 

H = f + V ext + V ee , 

where T is the operator of the kinetic energy, V ee the electron-electron inter- 
action, and l^xt the external potential. 

10.1.1 Hohenberg— Kohn Theorems 

Theorem 1 

The external potential V ext is, except for a trivial constant, uniquely defined 
by the electron density p{r). 

Proof: 

Assumption: There are two potentials V and V compatible with the electron 
density p{f). 

1) p — >V — >H — >ty 

2) p — >V' — > H' — >V' 
Using the variation principle for wavefunctions we get 

Eo<(V'\H\V') = ($'\H'\y') + (y'\H-H'\^) 
= E'o + J p{r)(y{r)-V\r))dr 
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and 

E' <(^\H'\^) = (*|#|#) + (V\H' -H\V) 

= E - J p(r) (V{r)-V\r))dr . 

We have used and as test functions with H' and H and have made 
use of the variation principle. After adding up the two equations we get the 
following contradiction 

Eq + E < E + E , 

and therefore our initial assumption must be wrong. Since we have used the 
variation principle this proof holds only for the ground state. 

Total energy 

E = E v (po) = T[ Po ] + J V cxt (r)p (r)dr + J[p ] + £ N cM 



Classical Coulomb-energy: 

Po(n)po(r 2 ) 



w = \fj 



\r 1 - r 2 \ 



dr\dr2 



Kinetic energy: T[p ] 

Non-classical electron-electron interaction energy: E N c [po] 
Only V cxt directly depends on the system, therefore we can write 

E = J V cxt (r)p (r)dr + F RK [p ] 

where F H k is a universal functional of the electron density: 

F RK [p]=T[p] + J[p]+E NC [p] . 



Theorem 2: 

For a trial density p(r), with the properties p(r) > and J p(r)dr = N 
(number of electrons) we always have E < Ey(p)- 

Proof: 

p ► V' —> H' —> 

HKl 
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We use 1^') as a trial function with the Hamilton operator H. We have 
(variation principal): 

E < = E v [p] 

and therefore 

E <E v [p] . 

This proof is based on the assumption that p — > V' is valid. But not all 

p(r) > with f p{r)dr = N have a corresponding potential. 

The set of all electron densities with a corresponding potential are called 

v-representable. 

The variation principle 

E = Mm^Evlp}} 
is only valid for v-representable electron densities. 

The proof can be expanded to iV-representable densities (Levi constraint 
search method) A iV-representable density is defined by: 

• calculated from an anti-symmetric wavefunction 

• fulfils 



p{r) > 0; J p{r)dr = N; J (Vp 1 



' 2 \ 2 dr 



< oo 



Plausibility arguments by Wilson 

1. Starting from p(r) 

2. Calculate the number of electrons N in the system 



N = J p{r)dr 

3. Find the position and charge of the nuclei (cusp-condition) 

d 



dr 



-p{r A ) 



= -2Z A p(Q) 

r A =0 



where p(r^) is the over all space directions averaged density. 
This defines the external potential V ext (r) = J2 A ^-r a \ 

4. N and V cxt define the Hamilton operator. 

5. With H we know \& and therfore all properties of the system can be 
calculated. 
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The electron density p(r) defines all properties of an iV-electron system. 
The Hohenberg-Kohn theorems give an existence proof for the functionals 
T[p] and £nc[p], as we U as a variation principle for the case of the known 
functionals. What is missing is the concrete form of the functionals. Many 
approximate functionals have been proposed. The accuracy, especially for 
the kinetic energy is for applications in chemistry not sufficient. In real 
application we will use another method. 

10.1.2 Kohn-Sham Theory 

We discuss a system with iV non- interacting electrons in a local potential V s : 

H s = f + V s . 



From the 1. Hohenberg-Kohn theorem we get 

E s [p]=T s [p] + J V s (r)p(r)dr 

and with the 2. Hohenberg-Kohn theorem it follows that 

5E s [p s ] = , 

This allows us to calculate the ground state density p s . 
The solution of the Schrodinger equation 



is a one-determinant wave function. This is due to the non-interacting elec- 
trons that allow for a simple product wavefunction for the solution. At the 
same time we also want the Pauli-principle to hold 

^ s = |$i(r)---$jv(r)| 

The one-particle density matrix for this system is 

N 

pl(r,r') = J2Mrmr') , 
i=i 

the density 

N 

i=i 
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and the kinetic energy 



T s [p\ = J **(r) (-^V 2 ) Mr)dr 



The functions <£>j are orthonormal. With this is the functional T s [p] for all 
densities of the form 



N 



p(r) = J2\Mr)\ 2 
i=i 

defined. 

Let's now turn to a system of N interacting electrons and an external poten- 
tial Vo- The corresponding energy functional is 

E Vo [p] = T[p] + J V (r)p(r)dr + J[p] + E NC [p] . 

Adding and subtracting T s [p] we get 

E Vo [p] = T s [p] + J V (r)p(r)dr + J[p] + T[p] - T s [p] + E NC [p] 

= T s [p] + / V {r)p(r)dr + J[p] + E xc [p] 



This defines the exchange-correlation functional E xc [p] 

E xc [p] = T[p] - T s [p] + E NC [p] ■ 

The corresponding potential is defined by the functional derivative 

5E xc [p] 
V xc [p]{r) = . 

5p{r) 

We assume that the optimal density po defined by 

5EvM\ P=Po = 

is non-interacting v-representable. Experience shows, that this is a reason- 
able assumption for physical systems. With that we have a unique potential 
V s with 

SEvAp]\ p=po =0 . 

Comparing the results we get 

V.(r) = V (r) + J y^dr' + V xc [ Po ](r) , 
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the definition of the Kohn-Sham potential. 

From the variation principle for the non-interacting system we get the eigen- 
value equation 



where the Lagrange multipliers q have been introduced for the orthogonality 
constraint of the Kohn-Sham orbitals $j. In analogy to the Hartree-Fock 
equation the potential V s depends on the solution of the Kohn-Sham equa- 
tions and the Kohn-Sham orbitals have to be calculated using an iterative 
procedure. 

Exactly like the Hartree-Fock orbitals the Kohn-Sham orbitals can be ex- 
panded in basis functions. We get a matrix representation of the Kohn-Sham 
operator 



where P is the one-particle density matrix of the Kohn-Sham orbitals. For- 
mally we can get the Kohn-Sham equations from the Hartree-Fock equations 
by replacing the matrix representation of the exchange operator by the ma- 
trix representation of the exchange-correlation potentials. 
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Lecture 11 



Density Functional Theory: 
Part II 



11.1 Exchange and correlation energy 
11.1.1 Density operators 

For a system with the wave function *ff N (xi, %2, • • • , ^at) (%i are combined 
spatial and spin coordinates) we define a density operator 7^ as 

7at(x' 1 , x' 2 ,..., x' n , X!, x 2 ,..., x N ) = ^ N (x[, x' 2 , x' N )^* N (x 1 ,x 2 , X N ) 

or in Dirac notation as a projection operator in position representation 

In = \^n)(^n\ ■ 

The trace of 7^ is one 

Tr( 7jv ) =/•••/ V N V* N dx! ■ ■ ■ dx N = 1 



and expectation values of operators A can be calculated from 

(A) = Tr( lN A) = Tr(A lN ) . 
For fermions we define reduced density operators of the form 

Specially we will need the two-particle density operator 

-f 2 (x 1 ,x 2 ,x 1 ,x 2 ) = / • • • / 7at dx 3 . . .dx 
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N 



'N 



and the one-particle density operator 

7i(x' 1 , xi) = N J ■■■ J 7^ dx 2 . . . dx N . 

Properties of 71 and 72 are 
. Tr( 72 ) = ^ 

• Tr( 7l ) = N 

• 71 and 72 are positive semi-definite and Hermitian and can be diago- 
nalised. 



7i(a; 1 ,xi)0i(xi) dxi = n^x^ 

j2(x 1 ,x' 2 ,x 1 ,X2)0i(x 1 ,x 2 ) dx 1 dx 2 = g$i{x\, x 2 ) 

4>i are the natural spin orbitals and Qi the natural spin geminals of the system. 
The eigenvalues n« and g^ are called occupation numbers of the orbitals and 
geminals, respectively. We have 

< rii,gi < 1 . 

and 

i 

N(N - 1) 

2> = — 2 — 

i 

Remember: 

Orbital one-electron function 
Geminal two-electron function 

In the basis of the eigenfunctions we can write the density operators as 

7i = 

i 

12 = J2^)&\ ■ 
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dsids 2 



Integration over spin variables: 

Pi(n,ri) = / 71(0:1, xi) , ds 1 
J s i= s i 

P2(r^2 ? r l> r 2) = j J 72(^1,4^1^2) , 1 = 

s 2 =s 2 

Pi and p 2 are spin-free density operators. Information about spin has been 
lost in p: 

Pi( r i> r i) = 7i(n<7,ri<7) 
Pi(K>n) = P?(n,ri) +Pi(ri,ri) 

Expectation values of operators that are not dependent on spin variables can 
be calculated by 



(01) = Tr(Oipi) = / Oi(ri)pi(r;,n) , dr[ 

J n=r 1 

(0 2 ) =Tv(0 2 p 2 ) = J J 2 (r 1: r 2 )p 2 (r' 1: r 2: r 1: r 2 ) 



ri=r 1 
r 2 =r 2 



dr\dr 2 



The electronic Hamilton operator is a two-particle operator, therefore we 
have 

E=(H} = Tr(Hp 2 ) = E(p 2 ) 

i.e. the total energy is a functional of the two-particle density operator. 
Electron density: 

p(r) = pi(r, r) 

p 2 (r,r') = p 2 (r,r ,r,r) 
11.1.2 Electron-electron interaction 

We write the diagonal part of the two-particle density operator in the form 



p 2 (n,r 2 ) = -p(ri)p(r 2 )[l + h(r u r 2 )] , 
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where h(ri,r 2 ) is the pair correlation function. The electron-electron inter- 
action energy is 



E Pf , — 



4p( r i)p( r 2)[l + M r i> r 2)] dr x dr 2 



\r 1 - r 2 \ 2 



= J ^ + 2 



Fi -r 2 | 
1 

\r 1 - r 2 \ 



- ] p{r 1 )p(r 2 )h(r 1 ,r 2 ) dr 1 dr 2 
p(ri)p xc (ri,r 2 ) dr x dr 2 . 



where p xc (ri,r 2 ) = p( r 2)h(r 1 ,r 2 ) is the exchange-correlation hole of an elec- 
tron at position r\. 
Sum rule: 



p(n) = 



N-l 

2 r l 

N-l J 2 



p 2 {r u r 2 )dr 2 
p(n)p(r 2 )[l + h(r 1 ,r 2 )}dr 2 



1 



TV- 1 

(N - l)p(n) = p( ri 



p(r 2 )dr 2 + J p{r 2 )h(ri,r 2 )dr 2 



N " + J p* c (ri,r 2 )dr 2 



From this we get 



J Pxc(ri,r 2 )dr 2 = -1 



11.1.3 One-determinant wave functions 



y S = \xi,X2,---,Xn) 



Reduced density operators 



-!(.;•,..(•.,.) : J y b (x 1: ...)V b *(x u ...)dx 2 ...dx N 

N 



l2{x'i,X 2 ,X 1 ,X 2 ) = ^ [71(^1,^1)71(4^2) - 7 1 (x' 1 ,X 2 )'Yi(x 2 ,X 1 ) 

72 is a simple function of 71. 
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Electron density: 

N 

i=i 

Pair correlation function: 

P2(r 1 ,r 2 ,r 1 ,r 2 ) = - [p(ri)p(r 2 ) - pi(ri, r 2 )pi(r 2 , n)] 
Be careful with spin! 
Pi(ri,r 2 )pi(r 2 , n) = p" (n, r 2 )p"(r 2 , n) + (n, r 2 )pf (r 2 , n) 

, / x Pi(n,7-2)pi(r 2 ,ri) 

/i ri,r 2 = / w n 

p(n)p(r 2 ) 

Sum rule: 

/ p{r 2 )h(r l ,r 2 )dr 2 = - [ p(r 2 ) ^^p}p^^ dr 2 
J J p(ri)p{r 2 ) 

Pi(r 1 ,r 2 )p 1 (r 2 ,r 1 )dr 2 



p(n) 

-7—- / 5^Xi(^i)xK r 2)xi(^)xl(n)^ 2 
^r^>i(n)x*(n) 



l 

p(n) 
-l 



The exchange-correlation hole of a Slater determinant fulfils the sum rule 
exactly. 

Invariance to unitary transformations: 



£/ CT is a unitary matrix (UU* = 1). Spin symmetry can not be broken by the 
transformation. 
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Example: 

i i 
ijk 
ijk 

3 

With pi is also P2 invariant and therefore the expectation values of the Hamil- 
ton operator and the total energy. 
Natural orbitals: 

7i = 



and 



7i 



and therefore rii = 1 and <fii = Xi- 
Idem-potency 



7i7i = J J 

= Z>i><&i 

= 7i 

11.1.4 Exchange and correlation energy in the Kohn— 
Sham method 

Kohn-Sham exchange and correlation functionals: 

E xc [p] = T[p] - T s [p] + V cc [p] - AP] 

We want to make a connection of Exc[p] with the pair correlation function 
from the electron-electron interaction energy. We consider Hamilton opera- 
tors where the operator V ee has been multiplied by a scaling factor A. For 
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A = 1 we get a system of fully interacting electrons and for A = a system 
of non-interacting electrons (similar to the Kohn-Sham reference system). 
The wave function to these operators will be denoted by ty x . From the 
Hohenberg-Kohn theorems we get for non-interactive v-representable p 

F x=1 [p]=F[p]=T[p] + V cc [p] 
F x=0 [p}=T s [p] . 

And therefore we have 

E xc [p]=T[p]-T s [p] + V cc [p]-J[p] 
= F x=1 [p] ~ F x=0 [p] - J[p] 




From a non-trivial derivation one gets finally 

dF x [p] 



and 



E XC [p]= /^A^lKel^)- J[p] 
Jo 

= j J ^7-/02 (n, r 2 ) dr x dr 2 - J[p] 

= \ l I — P( r i)p*c{ri,r 2 ) dr x dr 2 
2 J J r 12 

where the averaged exchange-correlation hole p xc is defined by 

PxcO"i,r 2 ) = p(r 2 )h(r 1 ,r 2 ) 
and the averaged pair correlation function h 
f 1 1 

J d\ p 2 (n,r 2 ) = -p(ri)p(r 2 )[l + h(r u r 2 )\ 
p 2 (ri,r 2 ) is the diagonal part of the two-particle density matrix calculated 



from the wave function ^/ x . 
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The averaged pair correlation function has the property 

h{r\,r 2 ) = h(r 2 ,ri) 

and the sum rule 

J Pxc(ri,r 2 )dr 2 = J p(r 2 )h(ri, r 2 )dr 2 = -1 . 

The exchange-correlation energy depends only on a spherically averaged 
Pxc(n,r 2 ) 

E xc[p] = \J drp{r) Ansdsp^(r,s) , 
where the spherically averaged exchange-correlation hole is defined by 

Px^( r > s) = J p xc (r, r')dr 
with the integration domain Q : \r — r'\ — s. The sum rule gets 

4tt j s 2 dsp s J(r,s) = -l . 

Definition of exchange and correlation energy in quantum chemistry: 

E xc [p]=Er[pn F ] + E^[p] 
E^ F : Exact energy - Hartree-Fock energy. 

E^ F '■ Hartree-Fock exchange energy calculated with Hartree-Fock orbitals. 
For the Kohn-Sham theory we will use a slightly different separation: 

E^[p]=E^[p]+E^[p} 

E^ s ■ Hartree-Fock exchange energy calculated with Kohn-Sham orbitals. 
The exchange hole for spin compensated systems can be calculated from 

/ x 1 |Pi(ri,r 2 )| 

and is 

J Px(ri,r 2 )dr 2 = -1 . 
The correlation hole is defined as the difference to p x 

p xc (n, r 2 ) = p x (ri,r 2 ) + p c (ri,r 2 ) 
and with that we have 



/ 



Pc(ri,r 2 )dr 2 = 



The exchange energy is identical to the Coulomb interaction of the electrons 
with a charge distribution with total charge one. The correlation energy is a 
Coulomb interaction of the electrons with a neutral charge distribution. 
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Lecture 12 



Density Functional Theory: 
Part III 



12.1 Local density approximation 

Exc[p] = J drp(r)e xc [p(r)] 

where exc( r ) = /(p( r )) is a function (not functional) of the density. 
Potential: 

Vxc{r) = exc(r) + p[r) 

dp{r) 

The approximation is to have at each point in space the energy density with 
the value of the homogenous electron gas of the local density. 
Separation of XC energy density 

exc(r) = e x(r) + e c (r) 
Hartree-Fock exchange for the homogenous electron gas (Dirac 1930) 

3 /3 ^ 1/3 



6x(p(r)) = -- ^-p 
Correlation energy for the homogenous electron gas: 

Quantum-Monte-Carlo calculations (Ceperley, Alder 1980, Ortiz, Ballone 
1994) 

Analytic function fitted to QMC values: 

r = [— 1 V3 
\_4np 
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1 + Piy/rl + (3 2 r s 
ec = A In r s + B + Cr s In r s + Dr s 



r s < 1 



7 C = -0.1423, /3i = 1.0529, (3 2 = 0.3334 
A = 0.0311, B = -0.048, C = 0.0020, D = -0.0116 

The formula for Exc m the local density approximation gives 



where /io(l r i — r 2|;p(^i)) is the averaged pair correlation function for the 
homogenous electron gas. Contrary to the exact formula we have here the 
density p(ri) as a factor and not p{r2)- The LDA exchange correlation density 
is 



but /?-o(| r i — r 2 1 ; p( r i) ) is n °t symmetric in r\ and r 2 . The exact function 
p xc is not spherical symmetric, but only the symmetric part is used in an 
exchange correlation energy calculation. Therefore, a compensation of errors 
can be expected from the LDA approximation. 

12.2 Gradient corrections 

Simple Taylor expansion don't work, but 'empirical' functionals of the type 



Pxc DA ( r i, r 2) = p(ri)/ioOi -r 2 \;p(r 1 )) 





is spherically symmetric and fulfils the sum rule 





give clearly improved results. 
Potential: 



SFxc 



3 



_d_ f £F XC 
dr a [5(V a p) 



Vxc{r) 



Sp{r) 



E 



a=l 
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12.2.1 Exchange functional of Becke (1988) 



E x = £^ DA [p(r)] - J F x [s]p 4/3 dr 
, IVpI 



p 4/3 



Fx[s]=(3 2- 1 l ? >- 



s 2 



1 + 6/9s sinh 1 s 

(3 optimised using Hartree-Fock exchange energies of noble gas atoms. 

P = 0.0042 

Typical error in the exchange energy Ex® A ~ 10% and £ x (Becke) w 0.1%. 

12.2.2 Correlation functionals of Perdew et al. 

(Perdew 1986, Perdew- Wang 1991, Perdew-Bourke-Ernzerhof (PBE) 1996) 
These functionals model the correlation hole based on exact properties with- 
out using empirical parameters. 
The PBE functional 



£ C PBE = J P{r) [e c + H(r s ,t)]dr 



1/3 / Q \ 1 / 3 

1 + At 2 



if(r s ,t) = 7 log l + -t 

7 



0,2 

with 



1 + At 2 + A 2 t 4 



^=-[exp(-e c / 7 )-l]- 1 

7 

7 = 0.031090690869, (5 = 0.066724550 

12.2.3 Lee- Yang-Parr (LYP) correlation functional (1988) 

This functional is based on an approximative formula for the correlation en- 
ergy of the Helium atom by Colle and Salvetti (1975). The basic assumption 
is 

P2(n,r 2 ) = Cp2 F (r 1} r 2 ) . 
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Then they replaced pJF with p 1 and after a partial integration one gets 



E 



LYP 



P 



1 + dx 



[1 + bC F e~ cx } 



dr 



+ / ab\Vp\ 2 x 



1 



1 + dx 24 



1 + 



cx + 



dx 



1 + dx 



where x = p x / 3 and 
C F = —(37T-)' ' << = 0.0 1918 



ll) (37r2)2/3 



0.132 



0.2533 



dr 



d = 0.349 



For |Vp| = (i.e. a homogenous system) the correct LDA value is not 
recovered. We have approximately 



E. 



LDA 



2£ c LYP (|Vp| = 0) 



The LYP functional has the special form 

F c (p,Vp) = |Vp| 2 G(p) 



12.3 Hybrid functionals 

Let's go back to the definition of the exchange-correlation energy in the 
Kohn-Sham method 

E xc [p] = T[p] - T s [p] + V cc [p] - J[p] 
= F x=l [p] - F x=0 [p] - J[p] 

and define the function E^ c as the non-classical Coulomb energy at the value 
A. With this we get 

E xc [p] = [ E x NC d\ . 
Jo 

For a value of A = (a non-interacting system) this energy is exactly known. 
It's the Hartree-Fock exchange energy, because the wavefunction for this 
system is a single Slater determinant. If we now approximate the end point 
of integration A = 1 with a common density functional, we can estimate the 
total integral with a weighted sum of the two end points. How to choose 
the weights is again an open question. There are many possibilities, some of 
them purely empirical and others based on theoretical considerations. 
The two most often used hybrid functionals are 
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Table 12.1: Dependence on basis set for some bond length [A]. LDA values 



Bond 


6-31G(d,p) 


6-311++G(d,p) 


Limit 


Experiment 


H 


-H 


-/0.765 


-/0.765 


-/0.765 


-/0.741 


H 3 C 


- CH 3 


1.513/1.105 


1.510/1.101 


1.508/1.100 


1.526/1.088 


H 2 C 


-CH 2 


1.330/1.098 


1.325/1.094 


1.323/1.093 


1.339/1.085 


HC 


-CH 


1.212/1.078 


1.203/1.073 


1.203/1.074 


1.203/1.061 



• B3LYP 

£f = (1 _ fl)jE LDA + aE HF + ^388 + ^LYP + (1 _ c)£? LDA 

where E^ 88 is the Becke exchange functional. The values for the pa- 
rameters are 

a = 0.20, 6 = 0.72, c = 0.81 

• PBE0 

^BE0 = E PBE + a25(E H F _ i ,PBE ) 

This is a parameter free functional. Neither the PBE functional nor 
the weight 0.25 are fit to any experimental data. 



12.4 Basis sets 

In density functional theory the same kind of basis sets can be used as in 
Hartree-Fock calculations. In principal the parameters would have to be 
newly optimised (e.g. for a 6-31G* type basis). In practice even individual 
parameters for each density functional would have to be generated. However, 
it turned out that such optimised parameters only have a minor effect on 
results and the expense needed for the optimisation is not justifiable. It is 
therefore common to use the same basis sets in density functional theory that 
are also used in wavefunction calculations. 

In Kohn-Sham density functional theory only simple orbitals have to be 
described. Therefore, the demand on the basis is comparable to other inde- 
pendent particle models like the Hartree-Fock model. Basis sets of double-^ 
quality with one set of polarisation functions are often sufficient to get good 
results. Triple-^ basis sets with two sets of polarisation functions are for 
most cases close to converged results. 
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Table 12.2: Averaged and maximal absolute error of the atomisation energy 
for a test set (G2) of molecules [kcal/mol]. Basis: 6-311+G(3df,2p) 



Functional 


RMS 


max error 


LDA 


36.4 


84 


PBE 


8.6 


26 


BLYP 


4.7 


15 


B3LYP 


2.4 


10 


PBEO 


3.5 


10 



12.5 Model calculations 

Comparing different functional for all kind of molecular properties is a favourite 
past time of all quantum chemists. There is a plethora of research papers on 
this topic. 

12.6 Coulomb energy 

Density functional programs are build exactly in the same way as Hartree- 
Fock programs. The only difference is that exchange integrals are replace 
by integrals with the potential from the XC functionals. Many DFT pro- 
grams differ in another point from Hartree-Fock programs. In Hartree-Fock 
program the Fock matrix elements for Coulomb and exchange are calculated 
together. 

Xa 

Computational time is dominated by the calculation of the two-electron in- 
tegrals. Therefore one uses 



F c <- 

/J.U 


- 4:P Xa (fiu, Xa) 




— AP^div, Xa) 




P»\(p>v, A<r) 


F x <- 


P v \(iw, Act) 




P/iaiH", Act) 




P va {liv, Xa) 
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to achieve a minimal number of integral calculations. In DFT calculation 
only the Coulomb part is needed, except for hybrid functionals. The Coulomb 
contribution has a special structure that can be used together with further 
approximations to increase the efficiency of the program. 
For the calculation of the Coulomb part of the Kohn-Sham matrix we need 
all integrals (/iv, Act) over all basis functions. This method scales like M 4 
where M is the number of basis functions. For large systems the number 
of integrals can be reduced drastically if only integrals with a significant 
contribution to the energy are calculated. All modern programs are using 
such screening procedures. It can be shown that asymptotically only M 2 
integrals have to be calculated. However, this procedure has a rather large 
pre-factor, i.e. the asymptotic scaling is only reached for very large systems. 
For this reason other approximate method for the calculation of the Coulomb 
energy are used. 

12.6.1 Resolution of identity method 

The electron density is defined as 



The product of two atom-centred basis functions is a new function 



The density is therefore described exactly in a basis of the size M 2 . 



The basic idea is now to replace this exact description of the density by 
an approximate description in a smaller basis, where the new basis is only 
slightly bigger than the basis used for the orbitals. 



p{r) = S ^ j P^X{nv)(j-) 

(pv) 




a 



With this we get for the contribution to the Kohn-Sham matrix 



a 





a 
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This method scales like M 3 and asymptotically again like M 2 , but it has a 
much smaller pre-factor. Special basis sets for the electron density have to be 
used that match the basis sets for the wavef unctions. Typically these basis 
sets have three times more functions than the wave function basis. 
We haven't yet discussed how the expansion coefficients c a are calculated. We 
can transform this problem into an optimisation. We are looking for those 
expansion coefficients which minimise the difference to the exact electron 
density with respect to a given metric. 

p( r ) = ^2c a x a {r) 

a 

A(r) = p{r) — p(r) 
n({c}) = J J A{r)0(r,r')A(r')drdr' 
c <- MIN {c} fi({c}) 
where 0(r,r') defines a distance metric. The solution of this problem is 

13 



with 



S a p = j J Xa (r)0(r,r') X p(r')drdr' 
T f3 = J J Xp{r)0{r,r')p{r')drdr' 



The operators 0(r,r') = 5(r — r'), 0(r,r') = \r — r'\ and 0(r,r') = 
have been investigated. Most of the programs use the Coulomb operator. 
These methods are called RI (resolution of identity) method, because they 
can be derived using the identity operator (for orthogonal basis sets). 

1 = J^HH 

a 

This method is available in many quantum chemistry codes z.B. Turbomole, 
Gaussian 2003, DeMON. 

12.6.2 Local expansion of density 

The local expansion of density are based on similar ideas as the RI methods, 
but instead of looking for a global approximation to the density separate 
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approximations to all pair densities are used. 

Pab{t) = p <xp®a{r)$p{r) 

af3 
a£A,f3eB 

Pab{t) « ^c^ r ) + ^ c vXv(r) 

fj,&A u&B 

The pair density of atoms A and B is expanded in a basis located at the 
atoms A and B. The advantage of this method is that the calculation of the 
expansion coefficients is done locally and is therefore independent of the size 
of the molecule. The disadvantage lies in the limitation of the expansion to 
the functions on the same atoms. To achieve good results large basis sets 
(with high angular momenta) have to be used. This method is used in the 
program ADF. 



12.7 Numerical integration 

With Gaussian functions as basis sets, all integrals needed in standard quan- 
tum chemistry calculations can be derived analytically. Integrals for the ex- 
change and correlation functionals needed in density functional theory cannot 
be calculated analytically even for special types of basis functions. They have 
to be calculated using numerical procedures. 

12.7.1 Voronoi polyhedra 

We investigate a general integral over all space 

I = J F{r)dr . 

In a first step we divide the integral into different parts using a weight func- 
tion wa(j-). The function wa{t) should be close to unity near nucleus A and 
close to zero near all other nuclei. The function should be continuous and 
differentiable and defined in all space. We require further that 

J2 w 4 r ) = 1 • 

A 
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With this we get 




A 



r 



= 5> 



A 

To calculate the weight functions we introduce confocal elliptic coordinates 
(A,/i, 4>). These coordinates are defined by 

n + r 2 



Dabei ist r\ die Distanz zum Zentrum 1 und r 2 die Distanz zum Zentrum 2, 
R12 der Abstand der beiden Zentren und (f> der Winkel um die Verbindungslinie 
der Zentren. Der Definitionsbereich der Koordinaten ist 



Die wichtige Koordinate ist /i. Fiir zwei Zentren A und B haben wir 

r A - r B 
Vab - —b 

tiAB 

und wir sehen, dass \iab den Wert —1 am Zentrum A und auf der Achse von 
A weg besitzt. \iab hat den Wert 1 am Zentrum B und auf der Achse von 
B weg. Definieren wir eine neue Funktion s{hab) 




A 



A 



[0,2vr] 
[l,oo] 
[-1,1] 




for — 1 < /jab < 
for < Hab < 1 



With this we can define the weigth functions 



Pa(t) = ] [ s(ii A b) 



B^A 
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The whole space is devided into Voronoi polyhedra with respect to the nu- 
clei. The discontinuity of the function s(/iab) at /iab = is not suited for 
numerical calculations. Therefore we choose 

where A m is taken such that s(— 1) = 1 and s(l) = 0. A typical value for m 
is 10. 

It is advantageous to use also the relative size of the atoms, i.e. we make a 
variable transformation such that different atoms get different weights. 

vab = I^ab + oab(1 — Ab) 

The parameter <iab is calcuated from the Bragg-Slater radii (Ra, Rb) of the 
atoms. 

uab 



o-ab = 



U AB- 1 



x-i 

uab = ——r 
X + 1 

_ Ra 

X ~ p 
Kb 

12.7.2 Radial integration 

We calculate the integrals I a using numerical integration, where the integra- 
tion is split into a radial part and an integration over the unit sphere. 

i 3 

There are many different integration methods for radial grids. The best meth- 
ods for this problem work in two steps. First the integral [0, oo] is mapped 
onto [—1,1] using a variable transformation. Possible transformations are 

a + x 

r = i] 



vT 



x 



2 



a + x 

r = rj- 

1 — x 



r] . a + 1 
In 



In 2 l-x ' 

where the transformations are normalised such that the midpoint of the x- 
interval is maped onto r — rj. The value of rj depends on the type of atom. 
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H 0.8 He 0.9 Li 1.8 Be 1.4 B 1.3 C 1.1 N 0.9 O 0.9 

F 0.9 Ne 0.9 

Na 1.4 Mg 1.3 Al 1.3 Si 1.2 P 1.1 S 1.0 CL 1.0 Ar 1.0 

K 1.5 Ca 1.4 Sc 1.3 Ti 1.2 V 1.2 Cr 1.2 Mn 1.2 Fe 1.2 

Co 1.2 Ni 1.1 

Cu 1.1 Zn 1.1 Ga 1.1 Ge 1.0 As 0.9 Se 0.9 Br 0.9 Kr 0.9 



Finally we use a Gauss quadrature for the integral on the interval [—1,1]. 

,1 N 

The form of the function W(x) depends on the orthogonal polynome that 
was used as a base for the integration method. 

Gauss-Legendre 
W(x) = 1, -1 < x < 1 

{j + l)P j+1 = (2j + l)xP j -jP j - 1 

Gauss-Chebyshev 

W(x) = -=L=, -Kx<l 
V 1 — x 2 

T j+ i = 2xTj - Tj-i 

12.7.3 Angular integration 

Two-dimensional integrals over the unit sphere 

/ / f(6,<l>)sin9dBd<j> 
Jo Jo 

can be efficiently performed using Lebedev grids. Lebedev grids with a cer- 
tain number of points allow for the exact integration of all spherical har- 
monics up to a certain value of L max . These grids are based on octahedral 
symmetry. The position of the special points and their integration weights 
have to be taken from a table. 
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Gitter 


^max 


Inred 


hot 


1 


3 


1 


6 


2 


5 


2 


14 


3 


7 


3 


26 


4 


9 


3 


38 


5 


11 


4 


50 


6 


15 


5 


86 


7 


17 


6 


110 


8 


19 


7 


146 


9 


23 


9 


194 


10 


29 


12 


302 


11 


35 


16 


434 


12 


41 


20 


590 


13 


47 


25 


770 


14 


53 


30 


974 
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Lecture 13 

Molecular Properties 



13.1 Wavefunction analysis 



Wavefunction ^ 



chemical concepts 
atomic charges 
oxidation states 
reactivity 
bonding 



13.1.1 Basisfunction based population analysis 

For \& = |$! • • • $7v| or 7(7-, r') = £V fi&i(r)$i(r') (single determinant wave- 
functions) we have 

i i 

Basis set expansion 

AO 

$ i( r ) = ^2c ai (j) a (r) 

a 

Inserting this into the density we get 

Pi(r) = fi\$i(r)\ 2 = fi^2c ai cpi(l) a (r)(l)p(r) 

a/3 
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Calulating the total number of electrons 

MO . MO . 

N = Pi{r) dr = J^/j^c^c^ J 4> a (r)4>p(r) dr 

MO 

= ^ ] fi ^ ] c ai c /3iS a /3 
i a/3 

a/3 

= Tr (PS) 
where we have used the density matrix 

MO 

i 

We will now look into schemes to partition the N electrons using (PS). 

Mulliken population analysis 

PaaS aa : fraction of electron in AO a AT „ 

„ „ , i , . iii ■ n JNow we sum over all a 

Paj3^aj3 '■ put | electron m a and ± electron m p 

associated with atom A. This only works with basis sets where the functions 

can be easily subdivided using atoms. 

AO AO 

PA = ^2 X/ P a fjSaf3 
aeA (3 

The total charge on atom A is then 

q A = Z A - p A 

Lowdin population analysis 

This population analysis is based on a slightly different calculation of the 
total number of electrons based on the invariance of the trace. 



N = Tr(PS) = Tr (PS^S^ 
= Tr (S^PS^ = Tr (P) 

= Vp 

/ 1 act 
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The number of electrons on atom A is 

AO 

Pa = ^2 ^ aa 

cteA 

The total charge on atom A is again 

Qa = Z A - p A 

These types of charges are 

• strongly depend on basis functions 

• can only be used for trends, not absolute values 

• use only with small basis sets (no diffuse functions) 

13.1.2 Electrostatic potential derived charges 

Die electrostatic potential of an electronic charge distribution p(r) and a 
collection of atomic cores with charge Za is 

W')=Erii--fri*' 

\r — Ra\ J \r — r\ 

The potential of a collection of atomic charges qA at the same positions Ra 
as the atomic cores is 

A 



\r-R A \ 



We now use the electrostatic potential to define the optimal atomic charges 
qA with 

• outside the charge distribution p{r) we require 

V ESP (r) = V£ SP (r) 

• charge conservation 

qA = total charge of molecule 

A 

• fix higher multipoles (dipole etc.) 

We achieve this using a least square fit with constraints 

Mh^ ^2 (^ESp(n) - ^Esp( r i)) 2 + constraints 

i 

These type of charges are often used in force fields and QM/MM calculations. 
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13.1.3 Density based charges 
Hirschfeld method 

Definition 



q A = Z A - / p(r) dr 

J Atom 



This requires a definition of an atomic volume. We start with a weighting 
function 

/ v PA{r) 

w A {r) = —j-^ 

where p A (r) is the electron density of an isolated atom and p(r) = J2 A f>A(r). 
The weighting function has the property 



J2w A (r) = 1 

A 

for any choice of r. Therefore we can write 

N — I p(r) dr — p(r) ^%(r) dr 
J J A 

= ^2 p( r ) w A(r) dr = ^p A 

A J A 

with the definition of p A 

Pa = J p(r)w A (r) dr 
Atoms in molecules (AIM) 

Bader analysis: based on the topology of p(r) (maxima, minima, saddle 
points). 

Maxima of p(r) are at nuclei = attractors of electronic charge. The gradi- 
ent of density at each point in space points to strongest attractor (basin of 
attraction = atom). 

13.2 Response properties 

Response of a molecule to external (internal) perturbation 
• external electric field (F) 
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Table 13.1: Example: Atomic charge of oxygen in water 



Basis set 


Mulliken 


Lowdin 


ESP 


AIM 


STO-3G 


-0.39 


-0.27 


-0.65 


-0 


89 


3-21G 


-0.74 


-0.46 


-0.90 


-0 


93 


6-31G** 


-0.67 


-0.44 


-0.81 


-1 


24 


6-311G(2d2p) 


-0.52 


0.00 


-0.74 


-1 


24 


6-311++G(2d2p) 


-0.47 


-0.12 


-0.76 


-1 


25 



• external magnetic field (B) 

• nuclear magnetic moment (I) 

• nuclear geometry (R) 

time-independent perturbation : static property 

time- dependent perturbation : dynamic property 

Total property: contribution from electrons (QM) and nuclei 

Taylor series expansion of energy 

Energy derivatives = "Properties" 
Example: External electric field 

E hlt = / p(r)V cxt (r) dr 



Assume electric field F = — is uniform over molecule (Multipole expan- 
sion) 

£int = qV - n ■ F - l -Q : F' - ■ ■ ■ 

q = net charge 
\i = dipole moment 
Q = quadrupole moment 
F' = field gradient 
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Field — > wavefunction — > dipole 

fi = permanent dipole 
a = polarisability 
P, 7 = hyper-polarisability 

The total energy is 

E(F) = | H - (j,F I *) 
with the total Hamiltonian 

H = H -fi F- aF 2 - ^PF 3 

Series expansion 

E(F) = E(0) H FH F H F H 

1 ' 9F , 2 OF 2 v 6 OF 3 

-HO -a 



General: 

Property oc ^ ^ 

n m Property 

Energy 

1 Electric dipole moment 
1 Energy gradient 

2 Polarisability 

2 Harmonic vibrational frequency 

1 1 Infrared absorption 

2 1 Raman intensity 



13.2.1 Derivative techniques 

Hamiltonian 

H = H + \P 1 + \ 2 P 2 

Expectation value 

E{\) = | H + APi + A 2 P 2 | *) 



138 



First derivative 
dE(\) 



< 



For real wavefunctions we have A — > then \1/ 



H | *) + (* | Pi + AP 2 I *> + (* I H 



9*, 
dA' 



<9P(A) 



<9A 



= (* I Pi I *o> + 2( 



A=0 



H | #o) 



A=0 



The wavefunction ^ depends on two different set of parameters, the wave- 
function expansion coefficients (c) and the basis functions (ip). The derivative 
wrt A is therefore 

_ dip dc 

For fully variational wavefunctions we have 

d^(a 1 ,a 2 ,...,a n 



da>i 

and we get the Hellmann-Feynman theorem 



d_ 



dH 



For variationally optimised wavefunctions (HF, DFT, MCSCF, but not MPx, 
CI, CCSD) we have 

= and _*0 



and therefore 

dE{\) 



dc 



(*o I I ^o) +2 ( 



A=0 



~d\ 



Hcllmann-Fcynman 



H | * > 



A=0 

Pulay term 



#0 

dA 
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Lecture 14 

NMR Chemical Shielding 



14.1 General 

Maxwell equations 

divE = 4np 
divB = 

1 OB 

TOtE = — — 

c at 

-j 4tt -> 1 <9£ 
rotfi = — J+ - — 

C COT 

Partial solution by introducing a vector field A with the property 

B = rot A 

A determines B uniquely, but A{B) is not unique. We can add a gradient 
field to A without change of the physical properties. 



A' = A + gradA 
rot A' = rot A + rot(gradA) = rot A 



A: gauge function 

Insert into Maxwell equation 



- ldB 

TOtE=— — 

c at 

E + -— = -grad$ 

C OT 
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$ : scalar potential with 

$' = $ + const. 

— * 

Special form of A: Coulomb gauge 

divA = 

Poisson equation 

A$ = -Anp 

. T 1 d 2 A 1 4tt 7 

-AA + -— + -grad— = —J 

& at 1 & at c 

In the non-relativistic limit (c — > oo) these equations simplify to 

A$ = -47rp 
Al= 

However, there is still a change in the gauge possible 

A' = A + gradA 

as we get for 

divA = div(A + gradA) 
= divgradA 
= AA 

This means that a change of gauge is possible for A that fulfill the Laplace 
equation 

AA = 

— * — * 

Special case of a homogeous external magnetic field B(r) = B 

A = ^B o x(r-R) 

= ^B x f- ^B x R 
1 

= -B x f + gradA 
X = ~(B xR)-f 

— * 

The position R is called the gauge origin. 
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14.2 Chemical shifts 

Nuclear spins couple to the local magnetic field B\ oc (Ra) 

• level split 

• transitions between levels can be induced by radiofrequency fields 

• measure energy difference oc B\ oc (Ra) 
Local magnetic field 

Bi oc (f) = B cxt (r) + B md (f) 
B cxt {f) = B 

Biot-Savart law 



1 f j(r) x (f — s) 

And s) = — / ^ ^ dr 
cj \r - s] 



where j(r) is the current density. We introduce the magnetic shielding con- 
stant 

B md (s) = -a(s)B 

and the chemical shift 

5 = Cstd - cr(s) 

— * 

The chemical shift is independent on the external magnetic filed B and the 
magnetic moment I of the nuclei. 

Calculation: For molecule and standard calculate a(s) at atomic positions, 
this needs the calculation of the current density j(r). 

Hamiltonian: 

H H 

p^p + A cxt 

— * 

For the Coulomb gauge and gauge origin R we have 

1 

2' 



H = H + A cxt p + ^A 2 cxt 



1 - 

A cxt = -B ext x (f-R) 



H(B cxt , R) = H + ^ ((f - R) x p) B cxt + - x (f - R)J 
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For one-determinant wavefunctions and closed shells we have 

# = |$i 

and for the current density 

N N 

j = i Y,(n^ k - (v*;)*o - 2A cxt n*k 
fc=i k=i 

14.3 Linear response for Hartree— Fock 

Fock operator 

h = h + ih 1 + C(B 2 cxt ) 

Hartree-Fock orbitals 

$ fc = $ fe +^ + O(5 c 2 xt ) 
The current density vanishes for B ext = 

We get for the first order current density 

N N 

k=l fc=l 

V 



V V 

paramagnetic part diamagnetic part 

j 1 is independent of the gauge origin R. The paramagnetic and diamagnetic 
parts depend on the gauge origin. This causes numerical problems if the two 
terms are not calculated with the same accuracy. This is called the gauge 
invariance problem of chemical shifts. It is effectively related to the problem 
of non-complete basis sets. The ground state and linear response orbitals 
have different needs for basis sets. 

Solution: distributed gauge methods 

IGLO individual gauge origins for every localised orbital 
GIAO gauge including atomic orbitals 
CSGT continous set of gauge transformations 
LORG localised orbital/local origin 



144 



Table 14.1: Example: chemical shifts for C calculated with different methods 
(CH 4 ) 



Basis 


common gauge origin 


IGLO 


GIAO 




C 


H 






TZ 


25.58 


62.01 


32.48 


32.57 


TZ+3p 


30.97 


41.22 


31.29 


31.35 


uncontracted 


30.98 


32.07 


31.20 


31.34 


extended 


31.23 


31.26 


31.24 


31.25 



14.3.1 Calculation of response orbitals 

We expand the linear response orbitals in the set of ground state orbitals. 

N 
9=1 

with Y qp = Y pq . Stationarity conditions for Hartree-Fock 

p _ q f i occupied 
m \o unoccupied 

Expansion of Fock matrix 
with 

FL = (®°i\h 1 \K + Y, H ^ Y ^ 

jb 

and 

H ia ,jb = $ijF° b - 5 ab F° + (aft, ij) - (aj, bi) 
The condition F} a = leads to a set of linear equations for Y of the form 

HY = h 

7i depends on <3>° and h depends on $° and h 1 . 

h l = -i\({?-R)*p)B cxt 

= ~((f-fl)xV)Bext 
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For canonical orbitals F pq = e p 5 pq and we get 

The induced current density is 

j 1 = 2 Yia (* °V$° - $° V<£>°) - (4xt x (r - £)) £ $°cj>° 

ia i 

Calculations of chemical shifts have been performed with Hartree-Fock, MP2, 
CCSD(T), and DFT electronic structure methods. For heavy elements rela- 
tivists corrections are needed. 
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